Accurate evolutions of inspiralling neutron-star binaries: prompt and delayed collapse to black hole 
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Binary neutron-star systems represent primary sources for the gravitational-wave detectors that are presently 
operating or are close to being operating at the target sensitivities. We present a systematic investigation in full 
general relativity of the dynamics and gravitational-wave emission from binary neutron stars which inspiral and 
merge, producing a black hole surrounded by a torus. Our results represent the state of the art from several points 
of view: (i) We use high-resolution shock-capturing methods for the solution of the hydrodynamics equations 
and high-order finite-differencing techniques for the solution of the Einstein equations; ( ii) We employ adaptive 
mesh-refinement techniques with "moving boxes" that provide high-resolution around the orbiting stars; ( Hi) We 
use as initial data accurate solutions of the Einstein equations for a system of binary neutron stars in irrotational 
quasi-circular orbits; (iv) We exploit the isolated-horizon formalism to measure the properties of the black holes 
produced in the merger; (v) Finally, we use two approaches, based either on gauge-invariant perturbations or 
on Weyl scalars, to calculate the gravitational waves emitted by the system. Within our idealized treatment 
of the matter, these techniques allow us to perform accurate evolutions on timescales never reported before 
(i.e. ~ 30 ms) and to provide the first complete description of the inspiral and merger of a neutron-star binary 
leading to the prompt or delayed formation of a black hole and to its ringdown. We consider either a polytropic 
equation of state or that of an ideal fluid and show that already with this idealized treatment a very interesting 
phenomenology can be described. In particular, we show that while higher-mass polytropic binaries lead to the 
prompt formation of a rapidly rotating black hole surrounded by a dense torus, lower-mass binaries give rise to a 
differentially rotating star, which undergoes large oscillations and emits large amounts of gravitational radiation. 
Eventually, also the hypermassive neutron star collapses to a rotating black hole surrounded by a torus. Finally, 
we also show that the use of a non-isentropic equation of state leads to significantly different evolutions, giving 
rise to a delayed collapse also with high-mass binaries, as well as to a more intense emission of gravitational 
waves and to a geometrically thicker torus. 

PACS numbers: 04.30.Db, 04.40.Dg, 04.70.Bw, 95.30.Lz, 97.60.Jd 



I. INTRODUCTION 

Little is required to justify the efforts in the study of bi- 
nary systems. Despite the simplicity of its formulation, the 
relativistic two-body problem is, in fact, one of the most chal- 
lenging problems in classical general relativity. Furthermore, 
binary systems of compact objects are considered one of the 
most important sources for gravitational-wave emission and 
are thought to be at the origin of some of the most violent 
events in the Universe. While some of the numerical diffi- 
culties involved in the simulations of such highly dynamical 
systems have been overcome in the case of binary black holes 
(BHs), numerical simulations of binary neutron stars (NSs) in 
general relativity have so far provided only rudimentary de- 
scriptions of the complex physics accompanying the inspiral 
and merger. Simulations of this type are the focus of this pa- 
per. 

Binary NSs are known to exist and for some of the systems 
in our own galaxy general-relativistic effects in the binary or- 
bit have been measured to high precision fl 0, . The in- 
spiral and merger of two NSs in binary orbit is the inevitable 
fate of close-binary evolution, whose main dissipation mech- 
anism is the emission of gravitational radiation. An important 
part of the interest in the study of coalescing systems of com- 
pact objects comes from the richness of general-relativistic ef- 
fects that accompany these processes and, most importantly, 



from the gravitational-wave emission. Detection of gravita- 
tional waves from NS binaries will indeed provide a wide vari- 
ety of physical information on the component stars, including 
their mass, spin, radius and equations of state (EOS) HI Ht]. 
Furthermore, NS binary systems are expected to produce sig- 
nals of amplitude large enough to be relevant for Earth-based 
gravitational-wave detectors and to be sufficiently frequent 
sources to be detectable over the timescale in which the de- 
tectors are operative. Recent improved extrapolations to the 
local group of the estimated galactic coalescence rates predict 
1 event per 3—10 years for the first-generation of interfer- 
ometric detectors and of 10 — 500 events per year, for the 
generation of advanced detectors |0] . 

There are three possible characteristic gravitational-wave 
frequencies related to the inspiral and merger of binary sys- 
tems. The first one is the frequency of the orbital motion of 
the stars in the last stages of the inspiral, before tidal dis- 
tortions become important. The second characteristic fre- 
quency is associated with the fundamental oscillation modes 
of the merged massive object formed after the onset of the 
merger. Numerical simulations in the frameworks of Newto- 
nian 0], post-Newtonian (PN) semi-relativistic |Ht] and 
fully general-relativistic gravity II lQfl have shown that, if a BH 
is not produced promptly, the frequency of the fundamen- 
tal oscillation modes of the merged object is between 2 and 
3 kHz, depending on the EOS and on the initial compactness 
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of the progenitor NSs. Finally, the third frequency is that of 
the quasi-normal modes (QNMs) of the BH, which is eventu- 
ally formed after the merger. 

The study of NS binary systems goes beyond the impact 
it has on gravitational-wave astronomy and is also finalized 
to the understanding of the origin of some type of 7-ray 
bursts (GRBs), whose short rise times suggest that their cen- 
tral sources have to be highly relativistic objects lUlll . Af- 
ter the observational confirmation that GRBs have a cosmo- 
logical origin, it has been estimated that the central sources 
powering these bursts must provide a large amount of energy 
(<~ 10 51 ergs) in a very short timescale, going from one mil- 
lisecond to one second (at least for a subclass of them, called 
"short" GRBs). It has been suggested that the merger of NS 
binaries could be a likely candidate for the powerful central 
source of a subclass of short GRBs. The typical scenario is 
based on the assumption that a system composed of a rotat- 
ing BH and a surrounding massive torus is formed after the 
merger. If the disc had a mass > 0.1 Mq, it could supply the 
large amount of energy by neutrino processes or by extracting 
the rotational energy of the BH. 

The understanding of GRBs is therefore an additional moti- 
vation to investigate the final fate of binaries after the merger. 
The total gravitational masses of the known galactic NS bi- 
nary systems are in a narrow range ~ 2.65 — 2.85 Mq and 
the present observational evidence indicates that the masses of 
the two stars are nearly equal. If this is the general situation, 
NSs in binary systems will not be tidally disrupted before the 
merger. As a result, the mass loss from the binary systems is 
expected to be small during the evolution and the mass of the 
merged object will be approximately equal to the initial mass 
of the binary system. Since the maximum allowed gravita- 
tional mass for spherical NSs is in the range ~ 1.5 — 2.3 Mq, 
depending on the EOS, the compact objects formed just after 
the merger of these binary systems are expected to collapse to 
a BH, either promptly after the merger or after a certain "de- 
lay". Indeed, if the merged object rotates differentially, the 
final collapse may be prevented on a timescale over which dis- 
sipative effects like viscosity, magnetic fields or gravitational- 
wave emission bring the star towards a configuration which 
is unstable to gravitational collapse. During this process, if 
the merged object has a sufficiently high ratio of rotational 
energy to the gravitational binding energy, it could also be- 
come dynamically unstable to nonlinear instabilities, such as 
the bar-mode instability lfl2l [l3ll . It is quite clear, therefore, 
that while the asymptotic end state of a binary NS system is a 
rotating BH, the properties of the intermediate product of the 
merger are still pretty much an open question, depending not 
only on the nuclear EOS for high-density neutron matter, but 
also on the rotational profile of the merged object and on the 
physical processes through which the object can lose angular 
momentum and energy. 

Several different approaches have been developed over the 
years to tackle the binary NS problem. One of these ap- 
proaches attempts to estimate the properties of the binary evo- 
lution by considering sequences of quasi-equilibrium config- 
urations, that is by neglecting both gravitational waves and 
wave-induced deviations from a circular orbit; this is expected 



to be a ve ry g ood approximation if the stars are well sera- 
rated 01 MM 

Other approaches have tried to simplify some aspects of the 
coalescence, by solving, for instance, the Newtonian or PN 
version of the h ydr odynamics equations (see (0, HI 
IH M El HI lH and references therein). At the same 
time, alternative treatments of the gravitational fields, such 
as the conformally-flat approximation, have been developed 
and coupled to the solution of the relativistic hydrodynamical 
equations ll34l l35tl . either in the fluid approximation or in its 
smooth-particle hydrodynamics (SPH) variant l36ll . Special 
attention has also been paid to the role played in these calcu- 
lations by the EOS and progress has been made recently with 
SPH calculations JlH. 

While all of the above-mentioned works have provided in- 
sight into the coalescence process and some of them represent 
the state of the art for their realistic treatment of the matter 
properties JHH], they are nevertheless only approximations to 
the full general-relativistic solution. The latter is however re- 
quired for quantitatively reliable coalescence waveforms and 
to determine those qualitative features of the final merger 
which can only result from strong-field effects. 

Several groups have launched efforts to solve the equations 
of relativistic hydrodynamics together with the Einstein equa- 
tions and to model the coalescence and merger of NS bina- 
ries The first successful simulations of bi- 
nary NS mergers were those of Shibata and Uryu j4lll42ll43ll . 
Later on, Shibata and Taniguchi have extended their numeri- 
cal studies to unequal-mass binaries providing a detailed and 
accurate discussion of simulations performed with realistic 
EOSs (see 
son et al. 

presenting results of binary NS evolutions using an adaptive- 
mesh-refinement (AMR) code. However, despite the high res- 
olution available with the use of AMR, no waveforms from the 
BH formation were reported in ref. l45ll over the timescales 
discussed for the evolutions. Finally, a set of simulations in- 
volving (among other compact binaries) also binary neutron 
stars have been recently presented by Yamamoto et al. using 
AMR techniques with moving boxes ll46tl . 

An aspect common to all the above-mentioned simulations 
is that, while they represent an enormous progress with re- 
spect to what was possible to calculate only a few years ago, 
they provide a description of the dynamics which is limited 
to a few ms after the merger. The work presented here aims 
at pushing this limit further and to provide a systematic in- 
vestigation of the inspiral, but also of the merger and of the 
(possibly) delayed collapse to a BH. Despite the fact that we 
do not account for the transport of radiation or neutrinos, our 
results benefit from a number of technical advantages (some 
of which are also shared by other groups): (i) The use of 
high-resolution shock-capturing methods for the solution of 
the hydrodynamics equations and high-order (i.e. fourth order 
in space and third order in time) finite-differencing techniques 
for the solution of the Einstein equations; ( ii) The use of adap- 
tive mesh-refinement techniques that provide higher resolu- 
tion around the orbiting stars; (Hi) The use of consistent initial 
data representing a system of binary NSs in irrotational quasi- 



44 ] and references therein). More recently, Ander- 
114511 have made an important technical progress by 
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circular orbits; (iv) The use of the isolated-horizon formalism 
to measure the properties of the BHs produced in the merger; 
(v) The use of two complementary approaches for the extrac- 
tion of the gravitational waves produced. Most importantly, 
however, our simulations can rely on unprecedented evolution 
timescales spanning more than 30 ms. 

Exploiting these features we provide, within an idealized 
treatment of the matter, the first complete description of the 
inspiral and merger of a NS binary leading to the prompt 
or delayed formation of a BH and to its ringdown. While 
our treatment of the matter is simplified with the use of ana- 
lytic EOSs, we show that this does not prevent us from repro- 
ducing some of the most salient aspects that a more realistic 
EOSs would yield. In particular, we show that an isentropic 
(i.e. poly tropic) EOS leads either to the prompt formation of 
a rapidly rotating BH surrounded by a dense torus if the bi- 
nary is sufficiently massive, or, if the binary is not very mas- 
sive, to a differentially rotating star, which undergoes oscil- 
lations, emitting large amounts of gravitational radiation and 
experiencing a delayed collapse to BH. In addition, we show 
that the use of non-isentropic (i.e. ideal-fluid) EOS inevitably 
leads to a further delay in the collapse to BH, as a result of the 
larger pressure support provided by the temperature increase 
via shocks. 

Our interest also goes to the small-scale hydrodynamics of 
the merger and to the possibility that dynamical instabilities 
develop. In particular, we show that, irrespective of the EOS 
used, coalescing irrotational NSs form a vortex sheet when the 
outer layers of the stars come into contact. This interface is 
Kelvin-Helmholtz unstable on all wavelengths (see, e.g. ll47ll 
and references therein) and, exploiting the use of AMR tech- 
niques, we provide a first quantitative description of this in- 
stability in general-relativistic simulations. 

Special attention in this work is obviously dedicated to 
the analysis of the waveforms produced and to their prop- 
erties for the different configurations. In particular, we find 
that the largest loss rates of energy and angular momentum 
via gravitational radiation develop at the time of the col- 
lapse to BH and during the first stages of the subsequent 
ringdown. Nevertheless, the configurations which emit the 
highest amount of energy and angular momentum are those 
with lower mass, because they do not collapse promptly to a 
BH, but instead produce a violently oscillating transient ob- 
ject, which produces copious gravitational radiation while re- 
arranging its angular-momentum distribution. We also show 
that although the gravitational-wave emission from NS bina- 
ries has spectral distributions with large powers at high fre- 
quencies (i.e. f > 1 kHz), a signal-to-noise ratio (SNR) as 
large as 3 can be estimated for a source at 10 Mpc if using the 
sensitivity of currently operating gravitational-wave interfer- 
ometric detectors. 

Many aspects of the simulations reported here could be im- 
proved and probably the most urgent among them is the in- 
clusion of magnetic fields. Recent calculations have shown 
that the corrections produced by strong magnetic fields could 
be large and are very likely to be present [see ref. ll48ll for 
Newtonian magnetohydrodynamical (MHD) simulations and 
ref. ll49ll50ll for a recent general-relativistic extension]. While 



we have already developed the numerical code that would al- 
low to perform the study of such binaries in the ideal-MHD 
limit 15111 . our analysis is here limited to unmagnetized NSs. 

The paper is organized as follows. In Section II we first 
summarize the formalism we adopt for the numerical solution 
of the Einstein and of the relativistic-hydrodynamics equa- 
tions; we then describe briefly the numerical methods we 
implemented in the Whisky code 152L 15311 . we outline our 
mesh-refined grid setup, and we finally describe the quasi- 
equilibrium initial data we use. In Sections III A and B we 
describe binaries evolved with the polytropic EOS and hav- 
ing a comparatively "high" or "low" mass, respectively. In 
Sections III C and D we instead discuss the dynamics of the 
same initial models when evolved with the ideal-fluid EOS, 
while Section III E is dedicated to our analysis of the Kelvin- 
Helmholtz instability. In Sections IV A and B we characterise 
the gravitational-wave emission for the case of the polytropic 
and ideal-fluid EOS, respectively. Finally in Sections IV C 
and D we report about the energy and angular momentum car- 
ried by the gravitational waves and their power spectra. In 
the Appendix, further comments on numerical and technical 
issues are discussed. 

We here use a spacelike signature (— , +, +, +) and a sys- 
tem of units in which c = G = M@ = 1 (unless explicitly 
shown otherwise for convenience). Greek indices are taken to 
run from to 3, Latin indices from 1 to 3 and we adopt the 
standard convention for the summation over repeated indices. 

II. MATHEMATICAL AND NUMERICAL SETUP 

A. Evolution system for the fields 

We evolve a conformal-traceless "3 + 1" formulation of 
the Einstein equations l54l l55l l56l l57tl . in which the space- 
time is decomposed into three-dimensional spacelike slices, 
described by a metric 7^, its embedding in the full spacetime, 
specified by the extrinsic curvature Kij, and the gauge func- 
tions a (lapse) and (3 l (shift) that specify a coordinate frame 
(see Sect. Ill Bl for details on how we treat gauges and lf58ll for 
a general description of the 3 + 1 split). The particular system 
which we evolve transforms the standard ADM variables as 
follows. The three-metric 7^ is conformally transformed via 

<t> = m det 7*j > 7ij = e'^Jij ( 1 ) 

and the conformal factor <f> is evolved as an independent vari- 
able, whereas jij is subject to the constraint dot 7^ = 1. The 
extrinsic curvature is subjected to the same conformal trans- 
formation and its trace triTy is evolved as an independent 
variable. That is, in place of A'y we evolve: 

K ee tr Kij = g ij K tj , Ay = (K^ - K) , (2) 

with tr Aij = 0. Finally, new evolution variables 

f 1 = f k t) k (3) 
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are introduced, defined in terms of the Christoffel symbols of 
the conformal three-metric. 

The Einstein equations specify a well-known set of evolu- 
tion equations for the listed variables and are given by 



(d t - Cp) 7y 



(dt - Cp) 4> 



-2aA 



--aK, 



(4) 



(5) 



discussion of the these points. 

Among the diagnostic quantities, we compute the angular 
momentum as a volume integral with the expression l60ll : 



j: 



vol 



Jjk 



1 



X^ lm jk A lm )e^d 3 x. (15) 



(d t - Cp) Aij = e-**[-DiD ja + a(R l3 - 8irSij)] TF 

+ a(KA ij -2A ik A k j ), (6) 
(d t - Cp) K = -D l D t a 



+ K 2 +An(p AnM + S) , (7) 



dtT 



2 



+ ^T l djl3 j - 2A ij dja + 2a(f l ]k A lk + 6A ij dj< 



(8) 



where Rij is the three-dimensional Ricci tensor, £)j the co- 
variant derivative associated with the three metric 7^, "TF" 



indicates the trace-free part of tensor objects and p A 
and Sij are the matter source terms defined as 

Padm = n a npT a(3 , 
Si = -ry ia npT a P, 



s, 



(9) 



where n a = (—a, 0,0,0) is the future-pointing four-vector 
orthonormal to the spacelike hypersurface and T a " is the 
stress-energy tensor for a perfect fluid (cf. eq. [T9V The Ein- 
stein equations also lead to a set of physical constraint equa- 
tions that are satisfied within each spacelike slice, 



K* - KaK» - 167rp ADM = 0, 



n = r {3) 1 

M l = Dj(K i: > - ^K) - 8ttS' 1 = 0, 



(10) 
(11) 



which are usually referred to as Hamiltonian and momentum 
constraints. Here = Rij"/ 1 -' is the Ricci scalar on a three- 
dimensional timeslice. Our specific choice of evolution vari- 
ables introduces five additional constraints, 



det 7y 
tvA 



1. 

13 — 0, 
p = jjkf 



(12) 
(13) 
(14) 



Our code actively enforces the algebraic constraints (fT~2T > 
and (fT3l . The remaining constraints, H, M. 1 , and ( fT4l . are not 
actively enforced and can be used as monitors of the accuracy 
of our numerical solution. See J59ll for a more comprehensive 



B. Gauges 

We specify the gauge in terms of the standard ADM lapse 
function, a, and shift vector, f3 l ll6lll . We evolve the lapse 
according to the "1 + log" slicing condition l62ll : 

d t a - {¥d t a = -2a(K - K ), (16) 

where K is the initial value of the trace of the extrinsic cur- 
vature and equals zero for the maximally sliced initial data 
we consider here. The shift is evolved using the hyperbolic 
T-driver condition 15911. 



dt(3 l - 0*8,? 



-aB l 



d t B l - frdjB 1 = d t T l - ftdjT 4 - ijB 1 



(17) 
(18) 



where 77 is a parameter which acts as a damping coefficient 
which we set to be constant (77 = 1). The advection terms on 
the rig ht-hand sides of these equations have been suggested 
in MM Ill- 
All the equations discussed above are solved using the 
CCATIE code, a three-dimensional finite-differencing code 
based on the Cactus Computational Toolkit [66]. A detailed 
presentation of the code and of its convergence properties have 
been recently presented in ref . J67ll . 



C. Evolution system for the matter 

An important feature of the Whisky code is the implemen- 
tation of a flux-conservative formulation of the hydrodynam- 
ics equations l68l l69l l70ll . in which the set of conservation 
equations for the stress-energy tensor T^ v and for the matter 
current density J* 4 , namely 



V„T' w = 0, VuJ" = 0, 



(19) 



is written in a hyperbolic, first-order and flux-conservative 
form of the type 



&q+3if W (q)=s(q), 



(20) 



where f W (q) and s(q) are the flux vectors and source terms, 
respectively lf7lll . Note that the right-hand side (the source 
terms) does not depend on derivatives of the stress-energy ten- 
sor. Furthermore, while the system ( |20b is not strictly hyper- 
bolic, strong hyperbolicity is recovered in a fiat spacetime, 
where s(q) = 0. 
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As shown by l69ll . in order to write system dT9b in the form 
of system d20"l l, the primitive hydrodynamical variables (i.e. 
the rest-mass density p, the pressure p measured in the rest- 
frame of the fluid, the fluid three-velocity v l measured by a 
local zero-angular momentum observer, the specific internal 
energy e and the Lorentz factor W) are mapped to the so- 
called conserved variables q = (D, S l , r) via the relations 

D = ^W Pl 

S i = v ^phW 2 v i , (21) 
t = ^( P hW 2 ~p)-D, 

where h = l+e + p/pis the specific enthalpy and 
W = (1 — 7y u 4 v-')~ 1 / 2 . Note that only five of the seven 
primitive variables are independent. 

In this approach, all variables q are represented on the nu- 
merical grid by cell-integral averages. The functions the q 
represent are then reconstructed within each cell, usually by 
piecewise polynomials, in a way that preserves conservation 
of the variables q IT72ll . This operation produces two values at 
each cell boundary, which are then used as initial data for the 
local Riemann problems, whose (approximate) solution gives 
the fluxes through the cell boundaries. A Method-of-Lines 
approach 17211 . which reduces the partial differential equa- 
tions d20b to a set of ordinary differential equations that can 
be evolved using standard numerical methods, such as Runge- 
Kutta or the iterative Cranck-Nicholson schemes IT731 l74ll . is 
used to update the equations in time (see ref. lf52ll for further 
details). The Whisky code implements several reconstruc- 
tion methods, such as Total-Variation-Diminishing (TVJD) 
methods, Essentially-Non-Oscillatory (ENO) methods 117511 
and the Piecewise Parabolic Method (PPM) [76]. Also, a va- 
riety of approximate Riemann solvers can be used, starting 
from the Harten-Lax-van Leer-Einfeldt (HLLE) solver 17711 . 
over to the Roe solver ll78ll and the Marquina flux formula 117911 
(see ll52l l53ll for a more detailed discussion). A comparison 
among different numerical methods in our binary-evolution 
simulations is reported in Appendix I A II 

Note that in order to close the system of equations for the 
hydrodynamics an EOS which relates the pressure to the rest- 
mass density and to the energy density must be specified. The 
code has been written to use any EOS, but all the tests so far 
have been performed using either an (isentropic) polytropic 
EOS 

p = Kp T , (22) 

e = P+Y&J' (23) 

or an "ideal-fluid" EOS 

p=(T-l)pe. (24) 

Here, e is the energy density in the rest frame of the fluid, K 
the polytropic constant (not to be confused with the trace of 
the extrinsic curvature defined earlier) and V the adiabatic ex- 
ponent. In the case of the polytropic EOS d22|, T = 1 + 1 /N, 
where N is the polytropic index and the evolution equation for 



t does not need to be solved. Note that polytropic EOSs ( f22b 
do not allow any transfer of kinetic energy to thermal energy, 
a process which occurs in physical shocks (shock heating). 
It is also useful to stress that by being isentropic, the poly- 
tropic EOS (1221 . far from being realistic, can be considered as 
unrealistic for describing the merger and the post-merger evo- 
lution. However, it is used here because it provides three im- 
portant advantages. Firstly, it provides one "extreme" of the 
possible binary evolution by being perfectly adiabatic. Sec- 
ondly, it allows us to use the same initial data but to evolve it 
with two different EOS. This is a unique possibility which is 
not offered by other (more realistic) EOSs. As we will com- 
ment below, it has allowed us to highlight subtle properties 
of the binary dynamics during the inspiral which were not re- 
ported before (see the discussion in Sec lIIICl ). Finally, by be- 
ing isentropic it provides the most realistic description of the 
inspiral phase, during which the neutron stars are expected to 
interact only between themselves and only gravitationally. 

In contrast to the polytropic EOS, when using the ideal- 
fluid EOS (l24l . non-isentropic changes can take place in the 
fluid and the evolution equation for r needs to be solved. 



D. Adaptive Mesh Refinement and Singularity Handling 

We use the Carpet code that implements a vertex- 
centered adaptive-mesh-refinement scheme adopting nested 
grids l80ll with a 2 : 1 refinement factor for successive grid 
levels. We center the highest resolution level around the peak 
in the rest- mass density of each star. This represents our rather 
basic form of AMR. 

The timestep on each grid is set by the Courant condition 
(expressed in terms of the speed of light) and so by the spatial 
grid resolution for that level; the typical Courant coefficient is 
set to be 0.35. The time evolution is carried out using third- 
order accurate Runge-Kutta integration steps. Boundary data 
for finer grids are calculated with spatial prolongation oper- 
ators employing third-order polynomials and with prolonga- 
tion in time employing second-order polynomials. The latter 
allows a significant memory saving, requiring only three time 
levels to be stored, with little loss of accuracy due to the long 
dynamical timescale relative to the typical grid timestep. 

In the results presented below we have used 6 levels of 
mesh refinement with the finest grid resolution of h = 
0.12 Mq = 0.177km and the wave-zone grid resolution of 
h = 3.84 Mq = 5.67 km. Our finest grid has then a res- 
olution of a factor of two higher compared to the one used 
in l44ll . where a uniform grid with h = 0.4 km was instead 
used. Each star is covered with two of the finest grids, so 
that the high-density regions of the stars are tracked with the 
highest resolution available. These "boxes" are then moved 
by tracking the position of the rest-mass density as the stars 
orbit and are merged when they overlap. In addition, a set 
of refined but fixed grids is set up at the center of the com- 
putational domain so as to capture the details of the Kelvin- 
Helmholtz instability (cf. Sect. llll Eb . The finest of these grids 
extends to r = 7.5 Mq = 11km. A single grid resolution 
covers then the region between r = 150 Mq = 221.5 km and 
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TABLE I: Properties of the initial data: proper separation between the centers of the stars d/M AriM ; baryon mass Mb of each star in solar 
masses; total ADM mass M ADM in solar masses, as measured on the finite-difference grid; total ADM mass M ADM in solar masses, as 
provided by the Meudon initial data; angular momentum J, as measured on the finite-difference grid; angular momentum J, as provided by 
the Meudon initial data; initial orbital angular velocity Qq ; mean coordinate equatorial radius of each star r e along the line connecting the two 
stars; ratio r e i /r e of the equatorial coordinate radius of a star in the direction orthogonal to the line connecting the two stars and r e ; ratio of the 
polar to the equatorial coordinate radius of each star r p /r e ; maximum rest-mass density of a star p m ax. The initial data for the evolutions with 
polytropic and ideal-fluid EOS are the same. Note that the asterisk in the model denomination will be replaced by "P" or by "IF" according 
to whether the binary is evolved using a polytropic or an ideal-fluid EOS. 
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r = 250 A/© = 369.2 km, in which our wave extraction is 
carried out. A reflection symmetry condition across the z = 
plane and a 7r-symmetry condition 1 across the x = plane are 
used. 

We have performed extensive tests to ensure that both the 
hierarchy of the refinement levels described above and the res- 
olutions used yield results that are numerically consistent al- 
though not always in a convergent regime. The initial data 
used for these tests refer to a binary evolved with an ideal- 
fluid EOS from an initial separation of 45 km and in which 
each star has a mass of 1.8 Mq\ such a mass is larger than the 
one used for the rest of our analysis (c/TableQ) and it has been 
employed because it leads to a prompter formation of a black 
hole, thus saving computational costs. 

The 2-norm of the typical violation of the Hamiltonian con- 
straint grows from < 10 -6 at the beginning of the simulations 
to < 10~ 4 at the end of the simulations. The convergence rate 
measured in the 2-norm of the violation of the Hamiltonian 
constraint is ~ 1.7 before the merger (i.e. the same conver- 
gence rate measured in the evolution of isolated stars l53ll ). 
but it then drops to ~ 1.2 during the merger. It is still unclear 
whether this difference is due to the generation of a turbulent 
regime at the merger (see the discussion in Sect. MI El and in 
ref. lf8llo or to a resolution which is close but not yet in a 
fully convergent regime. Tests showing the convergence rate, 
the conservation of the mass (baryonic and gravitational) and 
angular momentum, as well as the consistency in the gravita- 
tional waves have been validated by the referee but, for com- 
pactness, are not reported here. 

The apparent horizon (AH) formed during the simulation 
is located every few timesteps during the evolution l82ll . Ex- 
ploiting a technique we have first developed when perfor ming 
simulations of gravitational collapse to rotating BHs 11831 18411 
and that has now been widely adopted by other codes, we do 
not make use of the excision technique 18511 . Rather, we add a 
small amount of dissipation to the evolution equations for the 
metric and gauge variables only and rely on the singularity- 



Stated differently, we evolve only the region {x > 0, z > 0} applying 
a 180-degrees-rotational-symmetry boundary condition across the plane at 
X = 0. 



avoiding gauge ( fl6b to extend the simulations well past the 
formation of the AH (note that no dissipation is added to 
the evolution of matter variables). More specifically, we use 
an artificial dissipation of the Kreiss-Oliger type 18611 on the 
right-hand sides of the evolution equations for the spacetime 
variables and the gauge quantities. This is needed mostly be- 
cause all the field variables develop very steep gradients in 
the region inside the AH. Under these conditions, small high- 
frequency oscillations (either produced by finite-differencing 
errors or by small reflections across the refinement or outer 
boundaries) can easily be amplified, leave the region inside the 
AH and rapidly destroy the solution. In practice, for any time- 
evolved quantity u, the right-hand side of the corresponding 
evolution equation is modified with the introduction of a term 
of the type £ dis ,(«) = —eh 3 dfu, where h is the grid spacing 
and e = 0.075 and is kept constant in space (see ref. 084I1 for 
additional information and different prescriptions). 

E. Initial data 

As initial data for relativistic-star binary simulations we use 
the ones produced by the group working at the Observatoire 
de Paris-Meudon 1I20I, l25ll . These data, which we refer to also 
as the "Meudon data ", are obtained under the simplifying as- 
sumptions of quasi-equilibrium and of conformally-flat spatial 
metric. The initial data used in the simulations shown here 
were produced with the additional assumption of irrotational- 
ity of the fluid flow, i.e. the condition in which the spins of the 
stars and the orbital motion are not locked; instead, they are 
defined so as to have vanishing vorticity. Initial data obtained 
with the alternative assumption of rigid rotation were not used 
because, differently from what happens for binaries consisting 
of ordinary stars, relativistic-star binaries are not thought to 
achieve synchronisation (or corotation) in the timescale of the 
coalescence l87tl . The Meudon initial configurations are com- 
puted using a multi-domain spectral-method code, LORENE, 
which is publicly available 18811 . A specific routine is used to 
transform the solution from spherical coordinates to a Carte- 
sian grid of the desired dimensions and shape. 

The binaries used as initial-data configurations have been 
chosen so as to provide the variety of behaviours that we 
wanted to illustrate and some of their physical quantities are 
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reported in Table |U Furthermore, since it is the least com- 
putationally expensive, we have chosen model 1. 62-45-* as 
our fiducial initial configuration. For this binary the initial 
coordinate distance between stellar centres in terms of the ini- 
tial gravitational wavelength is d = 0.09 A GW , where A GW = 
7r/f2o is the gravitational wavelength for a Newtonian binary 
of orbital angular frequency Hq. For evolutions that employ 
a polytropic EOS, the polytropic exponent is T = 2 and the 
polytropic coefficient K = 123.6 = 1.798 x 10 5 g _1 cm 5 s~ 2 . 
(Ref. |01 has recently made the useful remark that a choice of 
r = 2.75 and K = 30000 leads to an EOS that fits well the su- 
pernuclear regime of the Shen-EOS at zero temperature l89ll . 
as well as yielding density profiles that are very similar to 
those obtained with that realistic EOS; unfortunately no initial 
data with this adiabatic exponent is available at the moment.) 



III. BINARY DYNAMICS 

In what follows we describe the matter dynamics of the bi- 
nary initial data discussed in the previous Section. To limit the 
discussion and highlight some of the most salient aspects we 
will consider two main classes of initial data, represented by 
models 1. 62-45-* and 1.46-45-*, respectively. These models 
differ only in the rest mass, the first one being composed of 
stars each having a mass of 1.625 M Q (which we refer to as 
the high-mass binaries), while the second one is composed of 
stars of mass 1.456 Mq (which we refer to as the low-mass 
binaries). 

Variations of these initial data will also be considered by 
changing, for instance, either the initial coordinate separa- 
tion (i.e. 60 km in place of 45 km) or the EOS (i.e. an ideal- 
fluid EOS or a polytropic one). Additional variations involv- 
ing, for instance, different mass ratios, will be presented else- 
where §01 



A. Polytropic EOS: high-mass binary 

We start by considering the evolution of the high-mass bi- 
nary evolved with the polytropic EOS, i.e. model 1.62-45-P in 
Table J] Fig.Q] in particular, collects some representative iso- 
density contours (i.e. contours of equal rest-mass density) on 
the (x, y) (equatorial) plane, with the time stamp being shown 
on the top of each panel and with the color-coding bar being 
shown on the right in units of g/cm 3 . 

The binary has an initial coordinate separation between the 
maxima in the rest-mass density of 45 km and, as we will dis- 
cuss more in detail later on, a certain amount of coordinate 
eccentricity and tidal coordinate deformation is introduced by 
the initial choice for the shift vector. As the stars inspiral, their 
orbital angular velocity increases and after about 2.2 orbits, or 
equivalently after about 5.3 ms from the beginning of the sim- 
ulation, they merge, producing an object which has a mass 
well above the maximum one for uniformly rotating stars, 
but which supports itself against gravitational collapse by a 
large differential rotation. Such an object is usually referred 



to as a hyper-massive neutron star or HMNS 2 . As the inspiral 
proceeds and the two NSs progressively approach each other, 
tidal waves produced by the tidal interaction become visible 
(cf. first and second rows of panels in Fig. [TJ and these are 
particularly large, i.e. of ~ 5%, for the high-mass binary and 
considerably smaller for the low-mass one (cf. Fig. [8}. 

This is shown in Fig. [2] which reports the evolution of the 
maximum rest-mass density normalized to its initial value. In- 
dicated with a dotted vertical line is the time at which the stars 
merge (which we define as the time at which the outer lay- 
ers of the stars enter in contact), while a vertical dashed line 
shows the time at which an AH is found. After this time the 
maximum rest-mass density is computed in a region outside 
the AH and therefore it refers to the density of the oscillat- 
ing torus. It is only a few orders of magnitude smaller. Note 
that before the merger the central rest-mass density not only 
oscillates but it also increases secularly, although at a much 
smaller rate (cf. also Fig.[6]i. 

As mentioned above, the merger takes place after about 
5 ms and the two NSs collide with a rather large impact pa- 
rameter. This reduces significantly the strength of the shocks 
which have been computed in the case of head-on colli- 
sions l92ll . but it also produces a considerable amount of shear, 
which could then lead to a series of interesting dynamical in- 
stabilities (see also the discussion in Sect. MI El l. Because of 
the adiabatic nature of the EOS, which prevents the forma- 
tion of physical shocks (i.e. discontinuities where entropy is 
increased locally) 3 , the HMNS produced at the merger is be- 
yond the stability limit for gravitational collapse, so that de- 
spite the high amount of angular momentum and the large de- 
gree of differential rotation 4 , it rapidly collapses to produce a 
rotating BH, at about 8 ms. 

More specifically, soon after the merger, the two massive 
and high-density cores of the NSs coalesce and during this 
rapid infall they experience a considerable decompression of 
~ 15% or more as shown in the small inset of Fig. [2] How- 
ever, after t ~ 6 ms, the maximum rest-mass density is 
seen to increase exponentially, a clear indication of the on- 
set of a quasi-radial dynamical instability, and this continues 
through the formation of an AH, which is first found at time 
t = 7.85 ms (see the last row of panels in Fig.[T]where the AH 
is shown with a thick dashed line, or Fig. [2] where the time of 



2 We recall that a HMNS is a star whose mass is larger than the maximum 
one allowed for a uniformly rotating model with the same EOS {i.e. the 
supramassive limit). For the T = 2 polytropes considered here this max- 
imum mass is M = 2.324 Mq (with an equivalent baryon mass of 
Mb = 2.559 Mq), while the maximum mass for a nonrotating model is 
M = 2.027 Mq (with an equivalent baryon mass of M b = 2.225 Mq). 
Clearly, all the models considered in Table Ulead to a HMNS after the 
merger. 

3 Note that very large gradients can nevertheless be produced because of the 
nonlinear nature of the hydrodynamic equations. These gradients, however, 
are not physical shocks in that no entropy is increased locally. 

4 Note that the HMNS is not axisymmetric and hence it is difficult to pro- 
vide a unique measure of the degree of differential rotation. On average, 
however, the angular velocity decreases of about one order of magnitude 
between the rotation axis and the surface. 
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FIG. 1: Isodensity contours on the (x,y) (equatorial) plane for the evolution of the high-mass binary with the polytropic EOS (i.e. model 
1.62-45-P in TableHJ. The time stamp in ms is shown on the top of each panel the color-coding bar is shown on the right in units of g/cm 3 
and the thick dashed line represents the AH. A high-resolution version of this figure can be found at l9lll . 
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FIG. 2: Evolution of the maximum rest-mass density normalized to 
its initial value for the high-mass binary. Indicated with a dotted ver- 
tical line is the time at which the stars merge, while a vertical dashed 
line shows the time at which an AH is first found and which is a few 
ms only after the merger in this case. After this time, the maximum 
rest-mass density is computed in a region outside the AH and there- 
fore it refers to the density of the oscillating torus. It is only a few 
orders of magnitude smaller. Note that the non-normalised value of 
the maximum rest-mass density at t = is 5.91 x 10 14 g/cm 3 (see 
Table HJ. The binary has been evolved using the polytropic EOS. 
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FIG. 3: Isodensity contours on the (x, z) plane highlighting the for- 
mation of a torus surrounding the central BH, whose AH is indicated 
with a thick dashed line. The data refers to the high-mass binary 
evolved with the polytropic EOS. A high-resolution version of this 
figure can be found at 19111 - 



appearance is marked by a dashed vertical line). 

This complex general behavior, namely the very small sec- 
ular increase in the central rest-mass density accompanied by 
small tidal oscillations, and the final decompression as the 
two NS cores merge into one, should help to clarify a long- 
standing debate on whether the NSs experience a compres- 
sion prior to the merger which leads them to c ollap se to a 
BH H [H |H or, rather, a decompression H3, HIH, as 
a result of the dynamical instability leading to the merger. 
Clearly, for the rather restricted set of stellar models which 
are close to the stability limit to BH collapse, the small sec- 
ular increase could lead to the formation of two BHs prior to 
the merger. 

After an AH is first found, the amount of matter outside 
of it is still quite large and, most importantly, it is the one 
with the largest amount of angular momentum. This leads to 
the formation of an accretion torus with an average density 
between 10 12 and 10 13 g/cm 3 , a vertical size of a few km 
but a horizontal one between 20 and 30 km (see evolution of 
Pmax in Fig.|2]after the AH). The torus has an initial rest mass 
of (M T )o — 0.04 M Q 5 , which however decreases rapidly to 



5 We define the initial mass of the torus as the rest mass outside the AH soon 
after the AH is first found. Note that such a measure could be ambiguous 
since the time of the first AH detection depends also on the frequency with 



become (M T )3 ms = 0.0117M Q only 3ms later. 

The dynamics of the torus are summarized in Fig. [3] which 
shows the isodensity contours on the (a;, z) plane; also in 
this case the time stamp is shown on the top of each panel, 
while the color-coding bar is shown on the right in units of 
g/cm 3 . Note that the panels refer to times between 13.2 ms 
and 16.7 ms and thus to a stage in the evolution which is be- 
tween the last two panels of Fig. Q] Other information on the 
properties of the merged object can be found in Table HI1 

Overall, the torus has a dominant m = (axisymmetric) 
structure but, because of its violent birth, it is very far from 
an equilibrium. As a result, it is subject to large oscillations, 
mostly in the radial direction, as it tries to compensate be- 
tween the excess angular momentum and the intense gravita- 
tional field produced by the BH. In doing so, it triggers quasi- 
periodic oscillations with a period of ~ 2 ms, during which 
the torus moves in towards the BH, accreting part of its mass. 



which the AH has been searched for and on the initial guess for the AH 
radius. To improve this notion and to give a measure that is indicatively 
comparable for different simulations, we take the values of the rest mass of 
the torus at the time at which the AH mean radius has reached the arbitrarily 
chosen value of 2.1. This mass should really be taken as an upper limit for 
the torus rest mass, since its value decreases considerably as the evolution 
proceeds and the torus accretes onto the BH. 
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FIG. 4: Evolution of the proper separation (top part) and of the co- 
ordinate separation (bottom part) for binaries with initial coordinate 
separation of either 45 or 60 km (i.e. models 1.62-45-P and 1.62-60- 
P in Table |IJ. Indicated with a dashed line is the proper separation 
for the binary starting at 45 km and suitably shifted in time. 



A behavior very similar to this one has been studi ed in de- 
tail in a number of related works JHH3, |H HI EM], in 
which the torus was treated as a test fluid. While the above- 
mentioned studies represent an idealization of the dynamics 
simulated here, they have highlighted that the harmonic dy- 
namics of the torus represent a generic response of the fluid to 
a quasi-radial oscillation with a frequency reminiscent of the 
epicy clic frequ ency for point-like particles in a gravitational 
field Il0llll02ll . Furthermore, because of the large quadrupole 
moment possessed by the torus and its large variations pro- 
duced by the oscillations, a non-negligible amount of gravi- 
tational radiation can be produced as a result of this process 
(see also the discussion of Fig.[T8ll. 

As mentioned in the Introduction, the existence of a mas- 
sive torus around the newly formed rotating BH is a key in- 
gredient in the modeling of short GRBs and the ability of 
reproducing this feature through a fully nonlinear simulation 
starting from consistent initial data is a measure of the matu- 
rity of simulations of this type. For compactness, we cannot 
present here a detailed study of the dynamics of the torus, of 
the variation of its mass and of the consequent accretion onto 
the BH. Such an analysis will be presented elsewhere ll90ll . but 
it is sufficient to remark here that the choice of suitable gauge 
conditions and the use of artificial viscosity for the field vari- 
ables allow for a stable evolution of the system well past BH 
formation and for all of the time we could afford computation- 
ally. 

Using the isolated-horizon formalis m Il03ll and its numer- 
ical implementation discussed in ref. Il04ll . we have mea- 
sured the final BH to have a mass M„„ = 2.99 M and spin 



7.3 Ml 



6.4 x 10 49 gcm 2 s _1 , 



thus with a dimen- 



sionless spin a = J BH /M| H = 0.82 (cf. Table HB. This is a 
rather surprising result when compared to the equivalent mea- 
sure made in the inspiral and merger of an equal-mass BH 
binary. In that case, it has been found that the final dimen- 
sionless spin is ag n ~ 0.68 for BHs that are initially non- 
spinning and increasing/decreasing for BHs that have spins 
para llel/anti-paralle l with the orbital angular momentum (see, 
e.g. lll05l Il06l Il07l Il08ll ). More specifically, the two initial 
BHs need to have a substantial spin, with a; n itiai — 0.45, 
in order to produce a final BH with afi na i ~ 0.82. On the 
other hand, the NSs have here little initial spin (they are es- 
sentially spherical besides the tidal deformation) and the lit- 
tle they have is anti-parallel to the orbital angular momentum 
(i.e. they counter-spin with respect to the orbital angular mo- 
mentum). Yet, they are able to produce a rapidly spinning 
BH. It is apparent therefore that the merger of two equal-mass 
NSs is considerably less efficient in losing the orbital angular 
momentum (or equivalently more efficient in transferring the 
orbital angular momentum to the final BH), thus producing a 
BH which is comparatively more rapidly spinning. 

An important validation of the accuracy of the simulations 
presented here can be appreciated when comparing the evo- 
lution of the same binary when evolved starting from differ- 
ent initial separations. More specifically, we have consid- 
ered high-mass binaries with initial coordinate separation of 
either 45 or 60 km (i.e. models 1.62-45-P and 1.62-60-P in 
Table|I]i and evolved them with a polytropic EOS. The results 
of this verification are summarized in Fig. |4j with the upper 
part reporting the evolution of the proper separation (continu- 
ous lines) and the lower one that of the coordinate separation 
(dotted lines) for binaries with initial coordinate separation 
of either 45 or 60 km (i.e. models 1.62-45-P and 1.62-60-P 
in Table J]). It should be remarked that the evolution of the 
latter binary is computationally much more challenging, with 
an inspiral phase that is about three times as long when com- 
pared with the small-separation binary. In particular, the stars 
merge at t ~ 18 ms, corresponding to ~ 5.5 orbits. This is 
to be compared with the ~ 2.2 orbits of 1.62-45-P and it is 
close to the limit of what is computationally feasible at these 
resolutions. 

The first thing to note in Fig. [4] is the remarkable differ- 
ence between the coordinate separation, which shows very 
large oscillations, and the proper separation, which instead 
shows only very little variations superposed to the secular de- 
crease. These are probably associated to a small but nonzero 
residual eccentricity like the one observed in binary BH sim- 
ulations i67ll . The oscillations in the coordinate separation, 
which have been reported also in ref. l45tl (cf. their Figs. 5 
and 7), are in our case clearly related to the gauge choice, as 
demonstrated by the evolution of the proper separation. This 
is also apparent when looking at Fig.|5j which shows the coor- 
dinate trajectory (dashed line) and the proper trajectory (solid 
line) of one of the two NSs in the high-mass binary start- 
ing from a coordinate separation of 45 km. While a certain 
amount of eccentricity is present also in the proper trajectory, 
this is rather small. 

A more careful analysis has revealed that the large oscil- 
lations in the coordinate separation are simply the result of 
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FIG. 5: Coordinate (dashed line) and proper (solid line) trajectory of 
one of the stars for the high-mass binary from a coordinate distance 
of 45 km (about 2.2 orbits). 



non-optimal gauge conditions. As mentioned in Sect. Ill El we 
import the initial data from the solution of the Meudon group, 
adopting the same shift vector (3 l computed for the quasi- 
circular solution. While this may seem like a reasonable thing 
to do, it actually introduces the oscillations commented above. 
We have also performed alternative simulations in which the 
shift vector has been set to be zero initially and then evolved 
with the gauge conditions ( 117V We have found that in this case 
also the coordinate separation is much better behaved and only 
very small oscillations are present (see also the discussion in 
Appendix IA 2\ . 

When concentrating on the evolution of the proper separa- 
tion it is clear that the binary starting at a large separation has 
a larger eccentricity, but also that most of it is lost by the time 
the stars merge. Indicated with a dashed line in Fig. |4]is also 
the evolution of the proper separation for the binary starting 
at 45 km when this is suitably shifted in time of ~ 13 ms; 
the very good overlap between the two curves is what one 
expects for a binary system that is simply translated in time 
and it gives a measure of our ability of accurately evolving bi- 
nary NSs for a large number of orbits. This provides us with 
sufficiently long waveforms to perform a first match with the 
PN expectations and also to establish the role played by the 
tidal interaction between the two NSs as they inspiral. Both 
of these studies will be presented elsewhere 1901. 

A comparison of the waveforms produced in these two 
simulations will be discussed in Sect. |IV A| but we show in 
Fig. [6] the evolution of the maximum rest-mass density nor- 
malized to its initial value for high-mass binaries with initial 
coordinate separation of either 45 or 60 km. For the large- 
separation binary, we observe a behavior very similar to that 
of the small-separation binary, as discussed for Fig. [2 namely 
the very small secular increase with superposed small tidal 
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FIG. 6: Evolution of the maximum rest-mass density normalized to 
its initial value for high-mass binaries with initial coordinate separa- 
tion of either 45 or 60 km. The vertical dashed lines denote the time 
at which an AH was found. The polytropic EOS was used for the 
evolutions. 



oscillations, the decompression as the two NS cores merge 
and the final exponential growth produced by the collapse to 
a BH. Note, however, that the two evolutions are not exactly 
the same and that small differences are appreciable both in the 
decompression phase and in the post-collapse phase which is 
dominated by the dynamics of the torus around the BH (Note 
that the different final values in /9 max are due to the different 
times at which the AH is first found and which do not coincide 
for the two runs; see also the footnote on page|9]l. 

Although a larger truncation error is to be expected in 
the case of the large-separation binary simply because of the 
larger integration time, we believe these differences are gen- 
uine and reflect the fact that the initial data used are not in- 
variant under time translation. Stated differently, the large- 
separation binary 1.62-60-P, when evolved down to a separa- 
tion of 45 km, will not coincide with the equilibrium solution 
1.62-45-P computed for a quasi-circular binary in equilibrium 
at 45 km. Because these differences are mostly in the internal 
structure, the deviations in the evolution become evident only 
at and after the merger and are essentially absent in the pre- 
merger evolution of both the central-density (cf. Fig. [6]i and of 
the waveforms (cf. Fig.|20]in Sect. lIV Al . 

Since considerations of this type have never been made be- 
fore in the literature and we are not aware of careful compara- 
tive studies of this type, our conclusions require further valida- 
tion. Work is now in progress to perform similar simulations 
with different polytropic indices. If the differences reported 
in Fig.|6]are indeed physical, they will also show variations in 
case stiffer or softer EOS are considered and they should in- 
deed disappear for perfectly incompressible stars. The results 
of this analysis will be reported in a future work l90ll . 
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FIG. 7: Isodensity contours on the (x, y) (equatorial) plane for the evolution of the low-mass binary with the polytropic EOS {i.e. model 
1.46-45-P in Table |IJ. The time stamp in ms is shown on the top of each panel, the color-coding bar is shown on the right in units of g/cm 3 
and the thick dashed line represents the AH. A high-resolution version of this figure can be found at f9lll . 
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FIG. 8: The same as in Fig. [2] but for the low-mass binary. Note 
that the merger time is essentially the same as for the high-mass bi- 
nary but there is a long delay in the collapse and the onset of quasi- 
harmonic oscillations in the HMNS. The binary has been evolved 
using the polytropic EOS. Note that the non-normalised value of the 
maximum rest-mass density at t — is 4.58 x 10 14 g/cm 3 (see 
Table©. 



0.5 

o 



0.5 



I 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

+ 3% " 


■ — 

i 


high-mass _ 


+ 3% " 


^ 

""""" — ~ — — 

J» * I 

tot 

_ J vol low-mass 

Jqws 





5 10 15 20 25 

t (ms) 



FIG. 9: Conservation of the total angular momentum for the high- 
mass binary (upper plot) and the low-mass one (lower plot). Indi- 
cated with different lines are the computed values of the volume- 
integrated angular momentum (solid line), of the angular momen- 
tum lost to gravitational waves (dotted line) and of their sum (dashed 
line). The dot-dashed line marks a 3% error. 



B. Polytropic EOS: low-mass binary 

We next consider the evolution of the low-mass binary 
evolved with the polytropic EOS, i.e. model 1.46-45-P in Ta- 
ble [H As for the high-mass binary, we first show in Fig. [7] the 
representative isodensity contours on the (x, y) plane, with the 
time stamp being shown on the top of each panel and with the 
color-coding bar being shown on the right in units of g/cm 3 . 
Note that because the evolution is different in this case, the 
times at which the isodensity contours are shown are different 
from those in Fig.Q] 

Since the mass difference with model 1 .62-45-P is less than 
10%, one expects that the orbital dynamics before the merger 
are essentially the same. Indeed this is what our simulations 
indicate and differences appear only as higher-order effects, 
such as in the strength of the tidal waves (see Fig. [7]). How- 
ever, despite the small difference in mass, the evolution after 
the merger is considerably different. This is nicely summa- 
rized in Fig.[8]which shows that the merger time is essentially 
the same as for the high-mass binary (i.e. ~ 5.3 ms), but the 
subsequent evolution does not lead to the prompt formation of 
a BH. Rather, the HMNS is still quite far from the instability 
threshold and undergoes a number of quasi-periodic oscilla- 
tions (cf. Fig. [8), which have almost constant amplitude in the 
central rest-mass density. 

A more careful analysis reveals that the core of the HMNS 
undergoes violent non-axisymmetric oscillations, with the de- 
velopment of an overall m = 2 deformation, i.e. a bar, as 
the system tries to reach a configuration which is energet- 



ically favourable through the rearrangement of the angular- 
momentum distribution. The bar-deformed star has an initial 
value of the ratio of the kinetic energy to the binding energy 
T/114 7 ! ~ 0.22, which remains roughly constant and only 
slightly decreases to 0.19 until the time of the collapse. As the 
bar rotates, it also loses large amounts of angular momentum 
through gravitational radiation and this is reported in Fig. [9] 
which shows the evolution of the angular momentum as nor- 
malized to the initial value. The top panel, in particular, refers 
to the high-mass binary, while the bottom one to the low-mass 
binary (a more detailed discussion of the losses of energy and 
angular momentum will also be presented in Sect. IIV Ob . In- 
dicated with different lines are the computed values of the 
volume-integrated angular momentum [solid line, computed 
with the integral (TT3T>1 . of the angular momentum lost to grav- 
itational waves (dotted line) and of their sum (dashed line). In 
both cases the slight secular increase is due to the truncation 
error and is at most of 3% over more than 20 ms (cf. dot- 
dashed line). A very similar figure can be made for the one of 
the ADM mass, whose conservation is even higher (the error 
is below 1%). For compactness we will not show such a figure 
here. 

Note that the loss of angular momentum is of ~ 5% of the 
total initial angular momentum during the inspiral and merger, 
but becomes much larger once the HMNS has been produced 
and while the bar-deformed core rotates. Indeed, in the case 
of the large-mass binary this loss increases to ~ 13% after 
the BH quasi-normal ringing, while it becomes as large as 
~ 22% for the low-mass binary. Overall, the post-merger 
evolution for the low-mass binary is rather long and spans 
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over ~ 16 ms. The inset in Fig. |8]shows that during this time 
the maximum rest-mass density oscillates but it also increases 
secularly of a factor of about 2. This is due to the fact that as 
the HMNS loses angular momentum, its centrifugal support 
is also decreased and thus it reaches more and more compact 
configurations. At one point the HMNS is past the threshold 
of the quasi-radial instability for the collapse to a BH, which 
takes place at ~ 20 ms (cf. Fig. [8), with an AH being found at 
t = 21.3 ms. 

Also in this case, a large amount of matter with sufficient 
angular momentum is found to be orbiting outside the BH 
in the form of an accretion disc. Differently from the high- 
mass binary, however, the torus here has a larger average 
rest-mass density (between 10 12 and 10 14 g/cm 3 ; see evo- 
lution of pmax in Fig. [8] after the AH), a larger extension in 
the equatorial plane (between 20 and 50 km) but a compa- 
rable vertical extension (below 10 km). It also has a larger 
baryon mass, which is initially (M T )o = O.IMq and be- 
comes (M T )3 ms = 0.0787 Mq after 3 ms (see footnote on 
page [9] and Table [II]). The dynamics of the torus are sum- 
marized in Fig. [TUJ which shows the isodensity contours on 
the (x, z) plane; note that the panels refer to times between 
21.4 ms and 27.4 ms and thus to a stage in the evolution which 
is between the last two panels of Fig. [7] A simple comparison 
between Figs.[3land[T0lis sufficient to capture the differences 
between the tori in the two cases and also to highlight that for 
a polytropic EOS the high-mass binary produces a lower-mass 
torus (cf. Table HT1 and see the discussion in Sect. lIVCl l. 

In analogy with what is seen for the high-mass binary, the 
torus has an overall axisymmetric structure and is far from 
equilibrium. As a result, it is subject to large oscillations, 
mostly in the radial direction, at a frequency close to the 
epicyclic one. A more detailed analysis of this will be pre- 
sented in a companion paper l90ll . 

Using again the isolated-horizon formalism we have es- 
timated that the final BH has in this case a mass M BH = 
2.60 Mq, spin J BH = 5.24M 2 = 4.61 x 10 49 gcm 2 s _1 and 
thus a dimensionless spin a = J BH /M 2 H = 0.76 (cf. TablelHI). 
Interestingly, the dimensionless spin is lower in the low-mass 
binary. 

It should also be remarked that the long time interval before 
which the collapse takes place has prevented previous studies 
from the complete calculation of the dynamics of NS binaries 
which would not lead to the p rompt formation of a BH. The in- 
vestigations of refs. l44j,l45ll49Tl . for instance, are limited to a 
few ms after the merger and should be contrasted with the evo- 
lutions reported here that cover a timescale of ~ 30 ms, also 
for the additional calculation of the gravitational waves. As a 
result, our simulations represent, within an idealized treatment 
of the matter, the first complete description of the inspiral and 
merger of a NS binary leading to the delayed formation of a 
BH. 



C. Ideal-fluid EOS: high-mass binary 

We now move on to discussing the dynamics of binary in- 
spiral and merger when the other EOS, the ideal-fluid one in 




-10 1 1 1 1 ' 

—40 -20 20 

X (km) 




-10 1 1 1 1 

-40 -20 20 

X (km) 




_ 10 l , 1 1 , L_ 

-40 -20 20 40 

X (km) 




-40 -20 20 40 

X (km) 

FIG. 10: Isodensity contours on the (x, z) plane, highlighting the 
formation of a torus surrounding the central BH, whose AH is indi- 
cated with a thick dashed line. The data refers to the low-mass binary 
evolved with the polytropic EOS (cf. Fig.PJ}. A high-resolution ver- 
sion of this figure can be found at 19111 . 



eq. $2M . is used. As discussed in Sect. Ill CI while this is an 
idealized and analytic EOS, it has the important property of 
being non-isentropic and thus of allowing for the change of the 
thermal part of the internal energy density (or, equivalently, of 
the temperature). As we will show in the remainder of this 
Section, this difference can lead to significant differences in 
the properties and dynamics of the HMNS produced by the 
merger. More specifically, we concentrate on the evolution of 
model 1.62-45-IF in Table [I] namely a binary in which each 
NS has a baryon mass of Mb = 1.625 Mq and an initial coor- 
dinate separation of 45 km. As for the previous binaries, we 
collect in Fig. QT|some representative isodensity contours on 
the equatorial plane. 

As one would expect from PN considerations (which sug- 
gest that finite-size e ffect s are expected at orders equal or 
smaller than the fifth lll09ll ). the bulk dynamics of the binary 
before the merger are essentially identical to the one already 
discussed for model 1.62-45-P and small differences are ap- 
preciable only in the low-density layers of the stars, where 
the different tidal fields cause comparatively larger amounts 
of matter to be stripped from the surface; this can be appre- 
ciated by comparing the second and third panels of Figs. Q] 
andQT] Indeed, this is a subtle point which is worth remark- 
ing: when using the ideal-fluid EOS, the evolution before the 
merger (i.e. during the inspiral) is not isentropic. This is be- 
cause small shocks are produced in the very low-density lay- 
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FIG. 11: Isodensity contours on the (x, y) (equatorial) plane for the evolution of the high-mass binary with the ideal-fluid EOS {i.e. model 
1.62-45-IF in Table|I]l. The time stamp in ms is shown on the top of each panel, the color-coding bar is shown on the right in units of g/cm 3 
and the thick dashed line represents the AH. A high-resolution version of this figure can be found at (9lll . 
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FIG. 12: Evolution of the maximum rest-mass density normalized to 
its initial value for the high-mass binary using the ideal-fluid EOS. 
Indicated with a dotted vertical line is the time at which the binary 
merges, while a vertical dashed line shows the time at which an AH is 
found. After this time, the maximum rest-mass density is computed 
in a region outside the AH. This figure should be compared with 
Fig.H 
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FIG. 13: Evolution of the coordinate volume-integral of the thermal 
part of the internal energy density. Note the secular increase in the 
thermal energy after the merger (vertical dotted line). The shaded 
area refers to a window in time in which the calculation of e t h via 
eq. d25t becomes inaccurate as a result of the steep gradients in the 
hydrodynamical variables developing inside the AH. 



ers of the stars as these orbit. These small shocks channel 
some of the orbital kinetic energy into internal energy, lead- 
ing to small ejections of matter (i.e. ~ 10 _6 Mq), and are thus 
responsible for the slight differences in the inspiral. We also 
note that these shocks would appear quite independently of 
the fact that the NSs are surrounded by an atmosphere as they 
represent the evolution of small sound waves that, propagating 
from the central regions of the stars, steepen as they move out- 
wards; we have checked that essentially identical results are 
obtained when changing the threshold for the atmosphere of 
one or more orders of magnitude (a discussion of this process 
for isolated stars evo lved within the Cowli ng approximation 
can be found in ref. ill lOll and in ref. Il 1 ill for the extension 
to a dynamical spacetime). 

Besides this small difference, the merger takes places at al- 
most the same time as for model 1.62-45-IF, namely after 
about 2.5 orbits, or equivalently after 5.8 ms from the begin- 
ning of the simulation. However, the post-merger evolution of 
the HMNS is considerably different. This is nicely summa- 
rized in Fig.Q~2l which reports the evolution of the maximum 
rest-mass density normalized to its initial value and which, af- 
ter the AH is found, refers to the region outside the AH. In 
this case shocks are allowed to form and the HMNS does not 
collapse promptly to a BH but, rather, undergoes very large 
oscillations with variations of 100% in the maximum of the 
rest-mass density (cf. Fig. [T2b , These oscillations are the re- 
sult of what appears to be a dynamical bar-mode instability 
which develops and is suppressed at least four times during 
the post-merger phase. More specifically, after the first ini- 



tial merger at t ~ 5 ms, the two stellar cores break up again to 
produce a bar-deformed structure, which rotates for more than 
a period before disappearing as the cores merge again. This 
process takes place four times and the merged object becomes 
increasingly more compact as it loses angular momentum and 
thus spins progressively faster. This behavior is clearly im- 
printed in the gravitational-wave signal as we will illustrate in 
Sect. UVBl 

Together with these large variations, the rest-mass den- 
sity also experiences a secular growth similar to the one al- 
ready discussed for the low-mass polytropic binary and, as 
discussed before, the increased compactness eventually leads, 
at t ~ 14 ms, to the collapse to a rotating BH. The use 
of the isolated-horizon formalism reveals that in this case 
the final BH has a mass M BH = 2.94 Mq, spin J BH = 
7.3 Mg = 6.4x 10 49 g cm 2 s _1 and thus a dimensionless spin 
a = J BH /M 2 H = 0.85 (cf. Table©. 

The explanation for this behavior in the post-merger phase 
and the appearance, also at high masses, of a delayed collapse 
to BH, can be found by considering the component of the spe- 
cific internal energy which is produced by the shock heating. 
This can be done by splitting the specific internal energy e into 
a cold component e co id = Kp /(r — 1) and into a thermal 
one e th , defined as 

p Kp r - 1 

eth = e_ecold = Kf— T)-T3T- (25) 

[Note that eq. (13) of ref. Il 12ll contains a typo for the expres- 
sion of e co id, which is corrected in expression (fZ5bl . Fig. [T31 
then reports the evolution of the coordinate volume-integral 
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of the thermal part of the internal energy density J v pethdx 3 , 
which increases secularly soon after the merger (vertical dot- 
ted line) [Indicated with the shaded area is the window in time 
in which the calculation of e t h via eq. d25T l becomes inaccurate 
as a result of the steep gradients in the hydrodynamical vari- 
ables developing inside the AH]. We recall that in the case of 
an isentropic EOS (such as the polytropic EOS), eth = and 
the quantity e/p r ~ 1 is constant and proportional to the spe- 
cific entropy of the system. However, with a non-isentropic 
EOS (such as the ideal-fluid EOS), entropy increases across 
shocks and the latter are clearly present during the merger, 
thus leading to a local and global increase of the specific in- 
ternal energy, namely of eth- As a result, soon after the merger 
(vertical dotted line in Fig.[T3l. the HMNS from an ideal-fluid 
high-mass binary can rely on an additional pressure support, 
which allows it to balance the gravitational forces at least for 
a few additional ms. Stated differently, the shocks produced 
at the merger are responsible for a local and global increase of 
the temperature, which will produce a global expansion of the 
HMNS and thus a reduction of its compactness. The overall 
smaller compactness caused by the increased internal energy 
can be appreciated by comparing the fourth and fifth panels of 
Figs.Q]and[n] 

A simple estimate for the temperature increase can be made 
by using the thermal part of the specific internal energy e t h 
and by neglecting the thermal energy due to radiation, so that 
eth = 3fcT/(2m„), where k is the Boltzmann constant and 
m n the rest mass of a nucleon. In this way the temperature is 
simply expressed as 



t=1 3.970 ms 



T = 



2m„ 



3fc(r - 1 

7.2174 x 10 12 
r - 1 



P 



P 



K. 



(26) 



Using ( 1261 1 it is then possible to estimate that the HMNS has 
an initial temperature of 5 x 10 10 K, which rapidly increases to 
5 x 10 11 K as the stellar cores merge. The additional shocks 
produced by the large oscillations in the post-merger phase 
can increase locally the temperature above these values, with 
maximum values that can reach 2 x 10 12 K. Clearly, at such 
large temperatures the radiative losses, either via photons or 
neutrinos, can become very important and lead to a qualita- 
tive change from the evolution described here. While first 
attempts of introducing the contribution of radiative losses 
in general-rela tivist i c ca lculations have recently been made 
(see, e.g. refs. JTT1 ITT1 ). we are still far from a mathemat- 
ically consistent and physically accurate treatment of these 
processes, which we will include in future works. For the 
time being it is sufficient to underline that, while it is clear 
that the inclusion of radiative processes will lead, quite gener- 
ically, to a decrease in the survival time of the HMNS after 
the merger, determining this time with any reasonable preci- 
sion will require not only the inclusion of radiative transport 
but also of a more realistic treatment of the EOS and of the 
scattering properties of the matter in the HMNS. 

In the absence of a more detailed calculation of the ra- 
diative losses, we can here resort to simpler back-of-the- 
envelope calculations to assess the importance of radiative 
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FIG. 14: Isodensity contours on the (x, z) plane highlighting the for- 
mation of a torus surrounding the central BH, whose AH is indicated 
with a thick dashed line. The data refer to the high-mass binary 
evolved with the ideal-fluid EOS (cf. Figs. [3] [H)]and the different 
vertical scales). A high-resolution version of this figure can be found 

at m. 



cooling in the post-merger phase when taking into account 
neutrino emission and diffusion. Let us therefore assume that 
the newly produced HMNS from a high-mass binary is ap- 
proximately spherical with an average radius of R HMNS ~ 
20 km, a mass of M HMNS ~ 3.2 M Q and thus an average 
rest-mass density which is essentially the nuclear rest-mass 
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density, i.e. p HMNS ~ p nuc ~ 3 x 10 14 g/cm 3 . We can 
now consider two different cooling processes acti ng ei ther 
via modified-URCA emission (see Chap. 1 1 of ref . llll5ln or 
through the more efficient direct- URCA emission ll 16ll . As- 
suming an initial average temperature of T HMNS ~ 10 11 K, 
the HMNS would cool down via modified-URCA processes 
to T HMNS - 10 10 (10 9 ) K in about 20 s (1 yr). On the other 
hand, if the cooling takes place through the much more ef- 
ficient direct-URCA processes, the cooling time would be 
~ 3 ms (1 min). Because the latter interval is smaller or com- 
parable with the ~ 9 ms elapsing in the present calculations 
between the formation of the HMNS and its collapse to a BH, 
we conclude that radiative losses in the HMNS would accel- 
erate its collapse to a BH only if direct-URCA processes take 
place. 

Quite predictably, also the merger of a high-mass binary 
evolved with the ideal-fluid EOS leads to the formation of 
a torus orbiting around the BH. With respect to the high- 
mass polytropic binary, however, the torus here has a differ- 
ent shape and a considerably larger vertical extension. Indeed 
the ratio of the vertical and horizontal sizes is ~ 0.5, while 
this was ~ 0.1 in the case of a polytropic EOS, irrespec- 
tive of the mass of the binary. Consequently, the measured 
initial rest mass of the torus is of a factor 6 larger than the 
one of the corresponding high-mass polytropic binary, namely 
(A/ T ) 3ms = 0.0726 M Q instead of (M T ) 3ms = O.O117M , 
3 ms after the first measure (cf. Table Hill. The average density, 
on the other hand, is considerably smaller (between 10 11 and 
10 12 g/cm 3 ). 

The dynamics of the torus are summarized in Fig. [14] which 
shows the isodensity contours on the (x, z) plane; note again 
that the panels refer to times between 14.0 ms and 22.4 ms 
and thus to a stage in the evolution which is between the last 
two panels of Fig. [14] A simple comparison between Figs. [3] 
[TOl and [T4l is sufficient to capture the differences among the 
tori in the three different cases considered so far. 

In view of the discussion made above on the increased in- 
ternal energy content produced by the shocks in the case of the 
ideal-fluid EOS, the formation of a vertically extended torus is 
not at all surprising, but the obvious response of the matter of 
the torus to a larger (thermal) pressure gradient in the vertical 
direction. Interestingly, the maximum rest-mass density of the 
torus does not show the typical harmonic behavior discussed 
so far in the case of the polytropic binaries and produced by 
the quasi-periodic oscillations in the radial direction. Rather, 
the maximum density shows a clear and monotonic decrease 
with time as a result of the accretion of the torus onto the BH 
(cf. Fig.[l~2]for £ > 14 ms). At the same time, the maximum of 
the internal energy in the torus is seen to increase (cf. Fig. [13] 
for t > 14 ms). Both the higher temperature and the geomet- 
rically thick shape of the torus produced in this case provide 
an important evidence that the merger of a massive NS binary 
could lead to the physical conditions behind the generation of 
a GRB. A more detailed analysis of the energetics and prop- 
erties of the torus (and in particular of its variability in time) 
is needed to further sup port this possibility and it will be pre- 
sented in a future work 119011 . 



D. Ideal-fluid EOS: low-mass binary 
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FIG. 15: Evolution of the maximum rest-mass density normalized to 
its initial value for the low-mass binary evolved using the ideal-fluid 
EOS. Indicated with a dotted vertical line is the time at which the 
binary merges. This figure should be compared with Fig. [2][8]and 
1121 of which maintains the same scale. 

Despite it being significantly different from the evolution 
of both the low-mass polytropic binary and of the high-mass 
ideal-fluid binary, the dynamics of the low-mass ideal-fluid 
binary is rather simple. In particular, the two NSs merge at es- 
sentially the same time as the corresponding high-mass ideal- 
fluid binary (i.e. t ~ 5.8 ms) and produce a HMNS which 
is however not sufficiently massive to collapse promptly to a 
BH. Rather, the HMNS undergoes a bar-mode instability pro- 
ducing an m = 2 deformation as the system tries to reach a 
configuration which is energetically favourable. 

Either as a result of the 7r-symmetry imposed (and which 
prevents the growth of the m = 1 mode) or simply because 
the HMNS is very close to the threshold of the bar-mode in- 
stability, the bar is seen to persist for the whole time the cal- 
culations were carried out, i.e. ~ 30 ms (see the discussion of 
ref. lfl3ll about under what conditions a bar-mode deformation 
is expected to survive over a longer timesc ale; recent addi- 
tional work on this can also be found in ref. ill 17ll ). Note that 
the bar deformation remains only approximately constant in 
time and that small oscillations in the central rest-mass den- 
sity can be measured. This is shown in Fig. [15] which reports 
the evolution of the maximum rest-mass density normalized 
to its initial value. Indicated with a dotted vertical line is the 
time at which the stars merge. This figure should be compared 
with Figs.[2l[8]and[T2l of which it maintains the same scale. 

During this rather long period of time (corresponding to 
~ 16 revolutions) the HMNS also loses large amounts of an- 
gular momentum through gravitational radiation (see discus- 
sion in Sects. lIVBl and llVCl ). As a result, the compactness of 
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TABLE II: Summary of the results of the simulations: proper separation between the centers of the stars d/M ADM ; baryon mass Mb of each 
star in solar masses; initial rest mass of the torus (M T )o (see footnote on page[5}; rest mass of the torus 3 ms after the appearance of the AH 
(M T )3 ms (actually 3ms after the time when the AH mean radius has reached the value 2.1, see footnote on page[9]l; mass of the BH M™, 
as computed in the isolated-horizon formalism; angular momentum of the BH J™ , as computed in the isolated-horizon formalism; BH spin 

parameter a = (J BH /M^ H ) IH , as computed in the isolated-horizon formalism; ratio of the ADM mass carried by the waves to the initial 
ADM mass; ratio of the angular momentum carried by the gravitational waves to the initial angular momentum. 
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FIG. 16: Left panel: Isodensity contours and velocity vector field (with the orbital component removed) on the (x, y) (equatorial) plane at a 
selected time soon after the merger. Note the presence of localized vortices in the shear layer between the two stars. Right panel: contours of 
the weighted vorticity p\ V x v\" (that is the rest-mass density multiplied by the module of the z component of the vorticity) for the same time 
shown in the left panel. This rendering highlights that in the shear layer the vorticity can be up to three orders of magnitude larger than in the 
bulk of the stars. Both panels refer to a high-mass binary evolved with the polytropic EOS. A high-resolution version of this figure can be 
found at (9l. 



the HMNS gradually increases and the central density shows 
the characteristic secular increase already discussed for the 
previous binaries (cf. inset of Fig. [TBI. The radiation-reaction 
timescale is in this case much longer: the HMNS is not very 
massive (its mass of ~ 2.6 M© is only sw 10% larger than the 
supramassive limit) and is more extended as a result of the 
increased internal temperature. As a result, the migration to 
the unstable branch and the collapse to a BH will occur much 
later than what was calculated and shown in Fig. [15] Using the 
latter to compute the growth rate of the maximum rest-mass 
density and assuming that the collapse to a BH is triggered 
when p max / p max (t = 0) — 2 (cf. Figs. l2l and [T21. we estimate 
that the collapse will take place at t ~ 110 ms. This timescale 
should be compared with the corresponding one (i.e. ^21 ms) 
obtained from the same initial data but evolved with a poly- 
tropic EOS. Clearly, the increase in the internal energy via 
shocks is responsible for this "long-delay" in the collapse to a 



BH. 

As a final comment we note that a timescale of ~ 110 ms 
is much longer than what is computationally feasible at the 
moment. As a result, the analysis of this binary will be lim- 
ited to a time interval of ~ 30 ms, which is, however, long 
enough to deduce its most interesting properties (see discus- 
sion in Sects. |WT3l and |T7Cl l. 



E. Vortex sheet and Kelvin-Helmholtz instability 

As mentioned above, when the two stars come into con- 
tact a vortex sheet (or shear interface) develops there where 
the tangential components of the velocity exhibit a discon- 
tinuity (i.e. the x and y components of the three-velocity in 
our setup). This condition is known to be unstable to very 
small perturbations and it can develop a Kelvin-Helmholtz in- 
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FIG. 17: Left panel: Maximum of the weighted vorticity p|V x v\ z on the equatorial plane normalized by the maximum of the rest-mass 
density p max during the evolution of the high (solid line) and low (dashed line) mass binaries evolved with the polytropic EOS. Indicated with 
a dotted vertical line is the time at which the binaries merge. Both the curves are plotted until the formation of an AH. Right panel: The same 
as in the left panel but for the ideal-fluid EOS. 



stabil ity, w hich will curl the interface forming a series of vor- 
tices 111 1811 . This is indeed what we observe in all our simu- 
lations, with features that are essentially not dependent on the 
mass or on the EOS used. 

In the left panel of Fig. [16] we show the isodensity con- 
tours and the velocity vector field on the equatorial plane for 
the high-mass binary evolved with a polytropic EOS at a time 
t = 6.091 ms when the presence of vortices is particularly 
evident. The density is shown in units of g/ cm 3 and in the 
bottom-right part of the plot an arrow is used as a reference 
for the values of the velocity. Furthermore, in order to high- 
light the formation of the shear interface, we have removed 
from the total velocity field the orbital angular velocity de- 
fined as the angular velocity of the stellar centers. The vector- 
field representation shows rather clearly that the vortex sheet 
goes from the bottom-left corner of the plot to the upper-right 
one. Along this sheet one can observe at least four main vor- 
tices, two of which are located at [x « ±7 km, y « ±5 km], 
while the other two are smaller and located at [x ss km, y as 
±2 km]. It is worth remarking that, because these smaller vor- 
tices have a scale of > 2 km ~ 1.3 M@, they are well captured 
by our resolution in the central regions which, we recall, is 
h = 0.12 Mq as 0.177 km. Because the employed numerical 
methods are very weakly dissipative on these scales, we be- 
lieve that our description of the Kelvin-Helmholtz instability 
is indeed accurate at the scales presented. Of course, differ- 
ent resolutions will either remove some of the vortices (as the 
resolution is decreased) or introduce new ones (as the resolu- 
tion is increased). In practice, we have found that a vortex of 
scale A is lost when the resolution used is h > 0.2A, proba- 
bly because the intrinsic numerical dissipation prevents their 
formation. 



A different and novel way of showing the presence of a 
vortex sheet and of the consequent development of a Kelvin- 
Helmholtz instability is offered in the right panel of Fig. [16] 
which shows the contours of the "weighted vorticity" on the 
equatorial plane, i.e. p|V x v\ z . Although this vector repre- 
sents the Newtonian limit of t he general-relativistic vorticity 
tensor ui^ = du,(huS) Qjj^] , it serves the purpose here of 
being proportional to the latter and also of simpler calcula- 
tion. Because the color-coding is made in a logarithmic scale, 
the right panel of Fig. [16] clearly highlights that the vorticity 
is not uniform in the merged object but that its value in the 
vortex sheet is up to three orders of magnitude larger than in 
the bulk of the stars. As stressed above, while both panels of 
Fig.[T6]refer to the high-mass polytropic binary, very similar 
results were obtained also for the low-mass binary or with the 
ideal-fluid EOS. 

To quantify the development of the Kelvin-Helmholtz insta- 
bility and measure its growth rate we have computed the maxi- 
mum of the weighted vorticity in the equatorial plane and plot- 
ted its time evolution in Fig.[l7j where it is also shown as di- 
vided by the maximum of the rest-mass density to remove the 
contribution due to the increase in p after the merger. Shown 
with different lines are the weighted vorticities for the high- 
mass binary (solid line) and for the low-mass binary (dashed 
line), evolved either with a polytropic EOS (left panel) or with 
an ideal-fluid EOS (right panel). Also indicated with a vertical 
dotted line is the time at which the two NSs merge, while the 
two vertical dashed lines refer to the times at which the AH is 
found in the two cases (no evolution is shown past this time 
as the measure of the vorticity becomes much more complex 
because of the turbulent motions in the torus). It is evident 
that after an initial growth of a factor of a few between t = 
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and t = 2 ms, probably produced by the transient away from 
the initial data, the weighted vorticity remains approximately 
small and constant. This stops at the time of the merger at 
t « 5 ms (cf. dotted vertical line) when the weighted vor- 
ticity grows exponentially of about two orders of magnitude. 
The Newtonian perturbative expectation for the growth rate is 
a ~ irv/A where v is the value of the velocity at the shear in- 
terface and A is the wavelength of the smallest growing mode; 
for v ~ 10~ 2 and A ~ 2 km, the measured growth rate is 
a ~ 10 3 s -1 and in reasonable agreement with the Newtonian 
expectation. 

The instability rapidly saturates when the two stellar cores 
merge; as a result, after ~ 2 ms from its initial development 
it reaches a quasi-stationary state. Note that the growth rate is 
essentially the same for the high- and low-mass binary and for 
the two EOSs (cf. the two panels Fig. [171) : however the evo- 
lution after the saturation is different for the different masses. 
The high-mass binaries collapse to a BH, while the HMNSs 
produced by the low-mass binaries hang on for a longer time, 
during which the instability persists at almost constant ampli- 
tude [for ~ 13 ms (cf. dashed line in Fig. ITTbl . 

As a final remark we note that, even if this instability 
is purely hydrodynamical, it can have strong consequences 
when studying the dynamics of binary NSs in the presence 
of magnetic fields. Indeed, as first shown by l48ll in New- 
tonian simulations and later briefly reported also by l49tl in 
general-relativistic calculations, in the presence of a magnetic 
field this instability leads to an exponential growth of the 
toroidal component even if the initial magnetic field is a purely 
poloidal one. In particular, it is reasonable to expect that even 
a moderate initial poloidal magnetic field of w 10 12 G can be 
increased up to values of order 10 15 G through this mecha- 
nism. Such high values of the magnetic fields are the ones 
presumed to be behind the phenomenology in magnetars, but 
are also thought to be the values necessary in order to extract 
sufficient energy from a system composed by a torus orbit- 
ing around a BH and power short hard GRBs. Work is now 
in progress for the investigation of this mechanism in fully 
general-relativistic MHD using the code presented in J5lll : re- 
sults of this investigation will soon be reported in a distinct 
article. 



IV. GRAVITATIONAL-WAVE EMISSION 

The accurate determination of the gravitational-radiation 
content of the simulated spacetimes represents a delicate and 
yet fundamental aspect of any modeling of sources of grav- 
itational waves; in view of this, we have implemented two 
different and equivalent methods to compute the gravitational 
waves produced by the inspiralling binaries. The possibil- 
ity of a comparison between the two methods and the cross- 
validation of the results provides us with additional confidence 
that the extracted waveforms are not only numerically accu- 
rate but also physically consistent. 

The first method uses the Newman-Penrose formalism, 
which provides a convenient representation for a number of 
radiation-related quantities as spin-weighted scalars. In par- 



ticular, the curvature scalar 

* 4 = -C a p 7 sn a mPn' y fri s 



(27) 



is defined as a particular component of the Weyl curvature ten- 
sor, C a p 1 s, projected onto a given null frame {I, n, m, m} 
and can be identified with the gravitational radiation if a suit- 
able frame is chosen at the extraction radius. In practice, we 
define an orthonormal basis in the three-space (r, 0,(f>), cen- 
tered on the Cartesian origin and oriented with poles along z. 
The normal to the slice defines a timelike vector t, from which 
we construct the null frame 
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We then calculate ^4 via a refor mulation of d27b in terms of 
ADM variables on the slice lll20ll. 



where 



Gij — F^ij KK^j 
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The gravitational-wave polarization amplitu des h + and h x 
are then related to ^4 by simple time integrals lll2lll 

h+ - ih x = * 4 , (31) 

where the double overdot stands for the second-order time 
derivative. 

The second and independent method is instead based on 
the measurements of the non-spherical gaug e-invariant per- 
turbations of a Schwarzschild BH (see refs. lfl22l [T2I \l2$\ 
for some applications of this method to Cartesian-coordinate 
grids). In practice, a set of "observers" is placed on 2-spheres 
of fixed Schwarzschild radius r s , derived from the coordinate 
(isotropic) radius via the standard formula 

- I • (32) 



2r i; 

where M = M ADM is assumed constant throughout the sim- 
ulation. On these 2-spheres we extract the gauge-invariant, 
odd-parity (or axial) current multipoles Q* m and even-parity 
(or polar) mass multipoles Q\ m of the metric perturba- 
tion Il25l Il26ll. Th e Q\ m and Q^ m variables are related to 
h+ and h x as Il27ll 



h + — ih y 
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Here _2Y tm are the s = — 2 spin-weighted spherical harmon- 
ics and (I, m) are the indices of the angular decomposition. 



A. Waveforms from polytropic binaries 

In what follows we illustrate and discuss the gravitational- 
wave signal produced by the inspiral and merger of the bina- 
ries discussed in Sect.|III]and we start by discussing the wave- 
forms produced by the binaries evolved with the polytropic 
EOS. 
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FIG. 18: Left panel: Retarded-time evolution of the real part of the I — m — 2 component of r^4 as extracted from a 2-sphere at a coordinate 
radius r = 200 Mq = 295 km for the high-mass binary. Indicated in the inset is the final part of the signal corresponding to the BH quasi- 
normal ringing. The merger takes place at (t — r) ~ 5.3 ms. Right panel: The same as in the left panel but shown in terms of the real part of 
the gauge-invariant quantity Q\ 2 - I n both cases the binaries have been evolved using the polytropic EOS. 
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FIG. 19: Comparison of the real part of the I = m = 2 component 
of r^4 for the high-mass binary evolved with the polytropic EOS 
when extracted at different radii: r = 160 Mq — 236 km (solid 
line), r = 200 M© = 295 km (dashed line), and r = 240 Mq = 
354 km (dotted line). 



Figure[T8] in particular, shows in the left panel the retarded- 
time evolution of the real part of the I = m = 2 component 
of r^4 as extracted from a 2-sphere at a coordinate radius 
r = 200 Mq = 295 km for the high-mass binary. Hereafter 
r = 200 Mq = 295 km will be the extraction radius for all 



the waveforms presented, unless specified differently. Indi- 
cated in the inset is the final part of the signal corresponding 
to the BH quasi-normal ringing. We recall that the merger 
takes place at (t — r) ~ 5.3 ms and that an AH is first found at 
(t — r) = 7.85 ms. The gravitational-wave signal during the 
inspiral is clearly very well captured and remarkably reminis- 
cent of the one observed in the man y binary B H simulations 
performed to date (see, for instance, 1 128lll29ll and references 
therein) and deviations from this type of waveforms are evi- 
dent only at (t—r) ~ 7 ms, when the HMNS starts its collapse 
to a BH. The ability of reproducing accurately the exponential 
decay of the quasi-normal ringing is often a good indication of 
having reached a sufficient level of accuracy as this involves 
the ability of measuring changes in the fields on the smallest 
possible physical scales (i.e. that of the horizon). The clean 
quasi-normal ringing shown in the inset shows that this is in- 
deed the case for the simulations reported here. 

It should also be added that, because the newly formed BH 
is not in vacuum but rather surrounded by a relativistic and 
accreting torus, the gravitational-wave signal should not be 
expected to be exponentially decaying to infinitesimal ampli- 
tudes during the ringdown. This explains the tiny but nonzero 
oscillations which can be seen after the ringdown and which 
are probably related to the accretion of matter onto the BH. 
A comparison with the results of re f. Ir%l WK llOOll or with 
the perturbative analysis of ref. [1 1 30ll could help to clarify the 
properties of this signal. 

The right panel of Fig. QjO on the other hand, shows the 
gravitational-wave signal in terms of the real part of the 
gauge-invariant quantity Q\ 2 - Because in this case the odd 
perturbations have zero real and imaginary part, the time evo- 
lution of the real (imaginary) part of Q\ 2 corresponds, modulo 
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FIG. 20: Comparison of the real part of the I = m = 2 compo- 
nent of r^i (upper panel) and of its amplitude (lower panel) for the 
high-mass binaries evolved with the polytropic EOS starting from 
an initial separation of 45 or 60 km. Indicated with a dashed line are 
the values after a time-shift. 



a constant coefficient, to the time evolution of the £ = m = 2 
component of h + (h x ). Note that the two waveforms are 
clearly different, but this is simply because they differ by two 
time derivatives [cf. eqs. ( 1311 and (13311. In fact, if a double 
time integral is made of the ^-waveforms, the corresponding 
curve lies on top of the one for Q l,. after a suitable normal- 
ization (see also Fig. 14 of ref. |67j where this was shown in 
the case of binary black holes). 

The comparison offered by Fig. [18] is useful to illustrate 
that, in contrast with what happens for binary BHs, the ampli- 
tude of the h + and h x polarizations does not increase mono- 
tonically in time but, rather, is reduced as the two NSs merge 
and as the HMNS collapses to a BH. Nevertheless, as we will 
comment in Sect. |IV CI the energy loss rate is largest during 
these stages (cf. right panel of Fig.[24l. 

Another important validation that the signal extracted cor- 
responds to gravitational radiation can be obtained by veri- 
fying that ^4 satisfies the expected "peeling" properties of 
the Weyl scalars, i.e. r 5 ~ n *f> n = const. This is illustrated in 
Fig. |T9l which compares the real part of the I = m = 2 com- 
ponent of when extracted at three considerably different 
radii: r = 160 Mq = 236 km (solid line), r = 200 M = 
295 km (dotted line), and r = 240 Mq = 354 km (dashed 
line) (the last radius is close to the outer boundary of our 
computational domain). Clearly, the overlap among the dif- 
ferent waveforms is very good both in phase and in amplitude 
and indicates that already at r ~ 150 Mq = 222 km gravi- 
tational waves can be extracted with confidence. (A similar 
figure can be built using the Schwarzschild perturbations and 
has not been shown here for compactness). 

It is interesting now to reconsider the impact that differ- 



ent initial separations of the same binary have on the emit- 
ted gravitational-wave signal. This aspect was already dis- 
cussed in Sect. IHI Al where the different dynamics were con- 
sidered, and nicely summarized in Figs. [4] and [6] We recall 
that the conclusions reached in Sect. IIII Al were that the differ- 
ences in the evolution of the large-separation binary 1.62-60- 
P and of its corresponding small-separation equivalent 1.62- 
45-P had to be found mostly in the internal structure and thus 
they were absent in the pre-merger evolution of both the cen- 
tral rest-mass density (cf. Fig. [6]) and the proper separation. 
A similar conclusion can be drawn also for the waveforms 
and we show in Fig. [20] a comparison in the real part of the 
I = m = 2 component of "J 4 (upper panel) for the high-mass 
binaries evolved starting from an initial separation of 45 or 
60 km. Note that the waveform for the 1.62-60-P binary con- 
tains more than 10 gravitational-wave cycles and is, therefore, 
the longest general-relativistic waveform computed to date. 

An equivalent view of this comparison is shown in the 
lower panel of Fig. [20] which reports instead the amplitude 
of x &4. Indicated with dashed lines in both panels are the val- 
ues after a suitable time shift. The good overlap in the inspiral 
phase is what is expected on PN grounds; however, a closer 
inspection also reveals that small differences do appear and 
these can then be used as a measure of the high-order PN cor- 
rections coming from compact binaries with finite size. More 
work and the use of long waveforms are necessary to study 
this further. 

We conclude this Section by discussing the gravitational- 
wave signal emitted by the low-mass binary and reported in 
Fig. [2T] Also in this case we show in the left panel the 
retarded-time evolution of the real part of the I = m = 2 
component of r\&4, while in the right panel the real part of the 
gauge-invariant quantity Q^- As mentioned in the previous 
Section, the HMNS has a prominent m = 2 bar deformation 
and gradually evolves towards a configuration which is un- 
stable to gravitational collapse through the emission of grav- 
itational waves. The loss of energy and angular momentum 
progressively reduces the centrifugal support and increases 
the compactness of the HMNS which, as a result, spins more 
rapidly. This is particularly clear in the evolution of ^4, which 
is shown in the left panel of Fig.[2T]and which exhibits the typ- 
ical increase in amplitude and frequency of the gravitational- 
wave signal. This runaway behavior ends at the time of the 
formation of the BH, which then rings down exponentially as 
shown in the two insets. A rapid comparison of Figs. [18] and 
Fig-Ellis sufficient to appreciate the marked differences intro- 
duced in the evolution of the binary by a different initial mass. 
In the following Section this comparison will be carried out 
also across different EOS s (cf. Fig.|23l. 



B. Waveforms from ideal-fluid binaries 

As mentioned when discussing the dynamics of ideal-fluid 
binaries, the significant differences that emerged both for the 
evolution of high- and low-mass binaries are reflected in their 
gravitational-wave emission. We recall that ideal-fluid bina- 
ries will experience a considerable increase of their internal 
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FIG. 21: Left panel: Retarded-time evolution of the real part of the I = m = 2 component of r^4 for the low-mass binary. Indicated in the 
inset is the final part of the signal corresponding to the BH quasi-normal ringing. The merger takes place at (t — r) ~ 5.3 ms. Right panel: 
The same as in the left panel but shown in terms of the real part of the gauge-invariant quantity Q\ 2 - I n bolh cases the binaries have been 
evolved using the polytropic EOS. 
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FIG. 22: Left panel: Retarded-time evolution of the real part of the I — m — 2 component of r^i for the high-mass binary evolved with the 
ideal-fluid EOS. Indicated in the inset is the final part of the signal corresponding to the BH quasi-normal ringing. Right panel: The same as 
in the left panel but for the low-mass binary. Here the inset does not refer to the BH quasinormal ringing; it highlights, instead, the periodic 
nature of the gravitational radiation after the merger. In both cases the merger takes place at (t — r) ~ 5.8 ms. 



energy (temperature) as a result of the shocks produced at the 
merger. As a result, a high-mass binary exhibits a delay in 
the collapse to BH of ~ 8 ms, which should be contrasted 
with the corresponding ~ 3 ms obtained for the same binary 
when evolved with a polytropic EOS. Similarly, a low-mass 
binary will show a much longer delay, which we estimated to 



be ~ 105 ms and which is to be contrasted with the corre- 
sponding ~ 16 ms obtained for the same binary when evolved 
with a polytropic EOS. 

This is nicely summarized in Fig. [22] whose left panel 
shows the retarded-time evolution of the real part of the 
£ = m = 2 component of r^4 for the high-mass binary. 
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FIG. 23: Left panel: Comparison in retarded-time evolution of the real part of the I = m — 2 component of r'1'4 for the high-mass binary 
when evolved with the polytropic or with the ideal-fluid EOS. Right panel: The same as in the left panel but for the low-mass binary. 



As commented in Sect. IIIICI the HMNS undergoes repeat- 
edly a dynamical bar-mode instability which develops and is 
suppressed at least four times during the post-merger phase, 
as the two stellar cores merge. The HMNS becomes increas- 
ingly more compact as it loses angular momentum and thus 
spins progressively faster. This behavior is clearly imprinted 
in the gravitational-wave signal and it is easy to distinguish 
the four stages of the bar development at times t ~ 8, 10, 12, 
and 14 ms, respectively. The last one is accompanied also by 
the gravitational collapse to BH and exhibits a well-captured 
quasi-normal ringing. 

The right panel of Fig. [22] on the other hand, refers to the 
low-mass binary and has a straightforward interpretation: the 
HMNS produced has a small m = 2 deformation and is still 
too far from the instability threshold to the collapse to a BH. 
Rather, the bar rapidly reaches an equilibrium configuration 
which persists over the 16 revolutions over which the calcu- 
lations were carried out. The resulting waveforms are pro- 
duced at twice the frequency of the revolution of the bar, i.e. at 
~ 2 kHz, and show a remarkably constant amplitude (cf. in- 
set in the right panel of Fig. |22] |. It is still unclear whether 
the stability of the deformation is the result of the bar being 
very close to the dynamical instability threshold or the result 
of the imposed 7r-symmetry, which prevents t he gr owth and 
coupling of the m = 1 and m = 2 modes lfl3l.ll 17ll . Clarify- 
ing this point will require calculations which are at least twice 
as expensive but it will be essential to determine whether the 
corresponding gravitational-wave spectrum will be character- 
ized by a large and predominant peak at ~ 2 kHz (cf. right 
panel of Fig. 1271). 

Figure [23] offers in its left panel a comparison in retarded- 



time evolution of the real part of the I = m = 2 component of 
r^4 for the high-mass binaries when evolved with the poly- 
tropic or with the ideal-fluid EOS (cf. left panels of Figs. [T8l 
andl22l. When shown in the same graph, it becomes much 
easier to appreciate the impact that the non-isentropic nature 
of the ideal-fluid EOS has on the dynamics of the merger 
and, most importantly, on the gravitational-wave emission. 
Clearly, when the waveforms from merging binary NSs will 
be detected, they will effectively provide the Rosetta stone 
for the deciphering of the stellar structure and EOS. In ad- 
dition, the comparison in Fig. [23] can also be used to gauge 
the possible range of behaviors that a more realistic treatment 
of the matter may yield. Both a polytropic and an ideal-fluid 
EOS, in fact, can be considered as the extremes of such a be- 
havior, with either a perfectly adiabatic evolution in which 
shocks (and hence shock heating) cannot occur, or with an 
evolution in which local increases of the temperature through 
shocks are allowed but cannot lead to radiative processes. Fur- 
thermore, because neutrinos or photons are expected to be 
trapped and would be able to leave the system only on dif- 
fusion timescales, any realistic EOS will lead to evolutions 
similar to those observed with the non-isentropic ideal-fluid 
EOS. 



Finally, the right panel of Fig. [23] is the same as in the left 
panel but for the low-mass binaries (cf. left panel of Figs. 1211 
and right panel of Fig. 122b . Also in this case the analogies 
and differences have a straightforward interpretation and un- 
derline the importance of considering the time between the 
merger and the collapse to BH as an important indicator of 
the properties of the binary. 
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C. Energy and Angular-Momentum Losses 

We have computed the energy and the angular momentum 
carried away by gravitational waves using the even and odd- 
parity perturbations, Q\ m an d Q* , respectively. The rate of 
energy loss, simply given by Il27ll 
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is shown in the right panel of Fig. [24] for all the low-mass 
and high-mass binaries considered here. In the left panel of 
the same figure we show the value of E GVJ normalized to the 
initial ADM mass of the system M AUM as a function of the 
retarded time t — r where r = 200 M@ = 295 km is the 
radius at which the waveforms were extracted. In both panels 
the solid line refers to the high-mass poly tropic model 1.62- 
45-P, the dashed line to the high-mass ideal-fluid case 1.62- 
45-IF, the dotted line to the low-mass polytropic binary 1.46- 
45-P, the dotted-dashed line to the low-mass ideal-fluid one 
1.46-45-IF and finally the long-dashed line to the high-mass 
polytropic model with an initial separation of 60 km, namely 
1.62-60-P. 

From the right panel of Fig. 124 lit is evident that all the mod- 
els have a first maximum in the energy emission rate soon after 
the merger. This initial increase in the emission rate is related 
to the last part of the inspiral phase, when the amplitude and 
the frequency of the gravitational-wave signal increase. Af- 
ter this first peak, however, the emission rate has a substantial 
drop, which is common to all the models and it is due to a 
very short (i.e. <C 1 ms) transition phase in which the devia- 
tions from axisymmetric are smaller. We now concentrate on 
describing the different dynamics of the different models after 
this initial common part, i.e. on the emission rate related to 
the evolution of the system after the merger. 

In the case of the two high-mass polytropic binaries, 
i.e. 1.62-45-P and 1.62-60-P, there is also a second peak in 
the energy emission at the time of the collapse and, except 
for the different times at which the merger and the subsequent 
collapse to BH take place, their profiles are very similar, with 
a total energy emitted which is ~ 0.01 M Aou . This second 
peak, which has an amplitude comparable to or higher than the 
first one, is simply related to the increase in amplitude and fre- 
quency of the gravitational waves emitted during the collapse 
(see also Fig.[20l>. In the case of the high-mass binary evolved 
with an ideal-fluid EOS, however, the emission rate exhibits 
four peaks after the merger and this is due to the different 
post-merger dynamics. As already discussed in Sect. lIIICl in- 
stead of collapsing promptly to a BH as the polytropic one, 
this system forms a bar-shaped HMNS with the high-density 
cores of the two NSs periodically merging and bouncing until 
sufficient angular momentum is carried away and the collapse 
starts. These periodic bounces and mergers of the two cores 
determine the several peaks seen in the emission rates. At 
the end, the total energy radiated through gravitational-waves 
is larger than the one emitted in the polytropic case and is 
~ 0.012 M ADM . 

For the two low-mass binaries, 1.46-45-P and 1.46-45-IF, 



on the other hand, the emission rate is always smaller than 
for the high-mass binaries, but it shows several peaks and for 
a longer time. This is related to the dynamics of the bar- 
deformed HMNSs that rotate for several stellar periods be- 
fore collapsing to BHs. As a result, even if the emission rate 
is smaller, the total energy emitted in gravitational waves is 
much larger and in the case of the low-mass polytropic bi- 
nary is ~ 0.018 M ADM at the time of the collapse, while 
for the low-mass ideal-fluid binary it can be estimated to be 
« 0.04 M ADM when extrapolating the time of the collapse to 
t w 110 ms (see discussion in Sect. HH St . 

The two panels in Fig.[24]are particularly useful to appreci- 
ate and quantify the differences that emerge among different 
binaries in the inspiral phase and, later on, in the post-merger 
phase. It is particularly instructive to consider the similarity 
in the evolutions of binaries having the same initial separation 
and mass, but different EOS, i.e. 1.62-45-P and 1.62-45-IF 
or 1.46-45-P and 1.46-45-IF. We recall that these sets of bi- 
naries have exactly the same initial data and hence the differ- 
ences during the inspiral are due uniquely to the role played 
by the EOS. As clearly shown in the left panel of Fig. [24] these 
differences are very small, so that 1.62-45-P and 1.62-45-IF 
have lost to gravitational waves essentially the same amount 
of mass at the time of the merger, although the latter actu- 
ally takes place at slightly different times (i.e. t ~ 5.3 ms for 
1.62-45-P and t ~ 5.8 ms for 1.62-45-IF; see the discussion 
in Sec. IIH Cl ). Because an identical comment also applies for 
1.46-45-P and 1.46-45-IF, we conclude that the EOS intro- 
duces major differences in the binary evolutions only after the 
merger. 

On the contrary, for binaries having the same EOS but dif- 
ferent masses (e.g. binaries 1.62-45-P and 1.46-45-P), also 
the evolution before the merger is different and can contribute 
to different post-merger evolutions (see the comment below 
on the angular-momentum losses). 

In a s imila r way, we have computed the angular-momentum 
loss as lfl27tl 



dt 



dt 



' ( III 



(QlmY 



(35) 



and, in analogy with Fig. [24] of which we use the same line- 
type convention, we show in the left panel of Fig.[25]the loss 
of angular momentum normalized to the initial angular mo- 
mentum of the system and in the right panel the loss rate. 

Overall, the angular- momentum losses and loss rates follow 
rather closely the behaviour already discussed for the energy, 
namely there is very little difference during the inspiral for bi- 
naries having the same mass. When looking more carefully, 
however, it is possible to note that the evolution of the angu- 
lar momentum shows small differences after about 3ms even 
for binaries with the same mass, e.g. 1.62-45-P and 1.62-45- 
IF. This time corresponds roughly to that of the first orbit, 
after which the non-isentropic evolution of model 1.62-45-IF 
will have changed the stellar structure only slightly but in a 
manner sufficient to produce a different emission of gravita- 
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FIG. 24: Left panel: Energy emitted in gravitational waves for the high-mass binary evolved with a polytropic EOS (solid line), for the low- 
mass binary evolved with a polytropic EOS (dotted line), for the high-mass binary evolved with the ideal-fluid EOS (dashed line) and for the 
low-mass binary evolved with the ideal-fluid EOS (dot-dashed line). Note that the largest amount of radiation comes from the low-mass binary 
whose emission has not been computed before. Indicated with a long-dashed line is the high-mass polytropic binary starting at 60 km. Right 
panel: The same as in the left panel but for the rate of the energy loss. Note that the largest burst of radiation is produced by the high-mass 
polytropic binary at the time of the prompt collapse to a BH. 
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FIG. 25: Left panel: The same as in Fig.|24]but for the orbital angular momentum normalized to its initial value (c^TablelUl- Right panel: The 
same as in Fig.[24]but for the rate of loss of orbital angular momentum. 



tional waves and hence a different loss of angular momen- 
tum. Indeed it is clear from the left panel of Fig. [25] that 
the merger begins when the two binaries have lost ~ 1% of 
their initial angular momentum but that this takes place at dif- 
ferent times for the two binaries, happening earlier for model 



1.62-45-P which is isentropic, more compact and with a larger 
quadrupole moment. 

More marked are the differences seen when comparing bi- 
naries differing only in the mass (e.g. binaries 1.62-45-P and 
1.46-45-P or binaries 1.62-45-IF and 1.46-45-IF). We recall 
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that these two sets of binaries essentially merge at the same 
time and it is then apparent from Fig. |25]that at the time of the 
merger the high-mass binary will have lost a larger relative 
amount of the initial orbital angular momentum. As a result, 
the matter orbiting outside the AH when this forms will also 
have a smaller amount of angular momentum and is therefore 
more likely to be more rapidly accreted. This explains why 
the high-mass polytropic binary 1.62-45-P produces a torus 
with a smaller rest mass than the low-mass polytropic binary 
1.46-45-P, both at the AH formation and after 3 ms {cf. Ta- 
bleHIli 6 . 

This behaviour indicates that, at least for binaries having the 
same EOS, the rate of loss of angular momentum during the 
inspiral phase plays an important role in determining the final 
mass of the torus and that the models that lose less angular 
momentum during the inspiral, hence comparatively low-mass 
binaries, are expected to have comparatively high-mass tori. 
This confirms what is already observed in ref. 114411 . 

Note, however, that such a simple conclusion is strictly true 
for binaries having the same EOS and when no radiative losses 
are taken into account. Under more generic conditions, how- 
ever, the EOS is also expected to play an important role and a 
representative example comes from comparing the high-mass 
binaries 1.62-45-P and 1.62-45-IF. In this case, in fact, the 
loss of angular momentum during the inspiral is essentially the 
same {cf. left panel of Fig. l25l l. but it is substantially different 
after the merger, with a loss of angular momentum which is 
at least 50% larger for the ideal-fluid binary. Yet, because of 
the increased pressure support the latter produces a torus with 
a mass which is ~ 7 times larger than the corresponding one 
for the polytropic binary. 



D. Gravitational-Wave Spectra and Signal-to-Noise Ratios 

We have also studied and compared the amplitudes and fre- 
quencies of the gravitational-wave signal produced by the dif- 
ferent models. In particular in Fig. [26] we plot the amplitude 
of the i = m = 1 component of the total gravitational-wave 

amplitude J h 2 + {t) + h 2 x (t) [where we neglect the contribu- 
tion of the spin-weighted spherical harmonic -2Y 22 in equa- 
tion d33lll. for four different binaries, all starting from an ini- 
tial separation of 45 km, as a function of the retarded time 
t — r, where r = 200 Mq = 295 km is the radius at which 
the signal was extracted. In particular, in the right panel we 
show the evolution of the gravitational- wave amplitude for the 
low-mass binaries 1. 46-45-* evolved using a polytropic EOS 
(solid line) and an ideal-fluid EOS (dashed line) while in the 
left panel we show the same but for the high-mass binaries 
1. 62-45-*. We recall that for all of these models the merger 
takes place after « 5 ms, which corresponds to the time when 
the amplitude reaches its maximum. The slight difference in 



6 Since we cannot follow the low-mass ideal-fluid binary till BH formation 
we cannot verify that this conclusion holds also for the ideal-fluid binaries, 
although we expect so. 



the position of these maxima between the polytropic and the 
ideal-fluid binaries is related to the difference in the time of 
the merger and is < 1 ms. 

Since the dynamics in the inspiral are very similar, the two 
high-mass binaries have a very similar and increasing ampli- 
tude, up to the merger. Note, however, that the increase is not 
monotonic and this is due mostly to the presence of a nonzero 
eccentricity. As commented in Sect. IIII Al a good part of the 
eccentricity is due to gauge effects (and is significantly re- 
duced when the shift vector is initially set to zero), but a small 
portion of it is also genuinely present in the initial data. Fortu- 
nately this spurious eccentricity has only a small impact in the 
power-spectral density (PSD) of the gravitational-wave signal 
and it is easy to isolate being it at ~ 4 times the orbital fre- 
quency. The evolution of the amplitude in the post-merger 
phase is rather different and, while it drops significantly for 
the polytropic binary, it remains at large values for the ideal- 
fluid binary as a result of the delayed collapse to BH; as we 
will comment later on, this will have an impact also on the 
detectability of this signal. 

The two low-mass binaries in the right panel of Fig. [26l also 
show a similar evolution up to the merger with an increase of 
the amplitude which is modulated by eccentricity and reaches 
its maximum at the merger. Of course, the maximum value 
reached in this case is lower than the one obtained in the high- 
mass cases. After the merger the amplitude is reduced by a 
factor of ~ 2 and remains to that level for the « 15ms which 
separate the merger and the collapse to a BH. In the case of the 
ideal-fluid binary, on the other hand, the post-merger ampli- 
tude is smaller and essentially constant for the whole time the 
simulation was carried out. As mentioned already, this binary 
is expected to collapse to a BH on a timescale of ~ 110 ms. 

We next consider the gravitational-wave emission in the 
frequency domain and for this we have computed the power- 
spectral density of the effective amplitude 

Mn = { hllf) f M , .36, 

where / is the gravitational-wave frequency and where 

poo 

h+,x(f)= / e 2 * ift h + , x (t)dt (37) 
Jo 

are the Fourier transforms of the gravitational-wave ampli- 
tudes /i+.x (t), built using only the largest I = m = 2 multi- 
pole. 

In Fig.|27]we compare the spectral distribution of the quan- 
tity h(f)f for the high-mass binaries (left panel) and the low- 
mass binaries (right panel) when evolved with the two EOSs. 
In both cases we use a solid line for the polytropic binaries and 
a dashed line for the binaries evolved with the ideal-fluid EOS. 
Also indicated in both panels with a vertical long-dashed line 
is the frequency corresponding to twice the initial orbital fre- 
quency /o = fio/(27r) where fo = 283 Hz for the low-mass 
binaries and fo = 295 Hz for the high-mass ones. These fre- 
quencies are representative of the signal at the beginning of 
the simulated inspiral and thus represent lower cut-off fre- 
quencies, below which the PSD is not meaningful. On the 
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FIG. 26: Left panel: Comparison of retarded-time evolution of the amplitude the I = m = 2 component of h = + ft^) f° r me high-mass 
binary when evolved with the polytropic (solid line) or with the ideal-fluid (dashed line) EOS; cf,, Fig. [23] left panel. Right panel: The same 
as in the left panel but for the low-mass binary; cf., Fig. [23] right panel. 




FIG. 27: Left panel: Comparison of the PSD of the i = m = 2 component of h(f)f for the high-mass binary when evolved with the 
polytropic (solid line) or with the ideal-fluid (dashed line) EOS; cf, Fig.[23] left panel. Right panel: The same as in the left panel but for the 
low-mass binary; cf, Fig. [23] right panel. Indicated with a vertical long-dashed line is twice the initial orbital frequency. 



other hand, the peaks in the PSDs observed at frequencies 
slightly larger than the orbital ones are very important as they 
refer to the power emitted during the inspiral. 

The PSD for the high-mass polytropic binary (left panel 
of Fig. |27| i is quite simple, as it shows, besides the inspiral 
peak, also a peak at / w 4 kHz, corresponding to the col- 



lapse of the HMNS (cf. left panel of Fig. [26]l. Note that the 
PSD shown does not include the frequency of the fundamen- 
tal QNM of the newly produced BH. Using the approximate 
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TABLE III: Signal-to-noise ratio (SNR) computed for different detectors assuming a source at 10 Mpc. The different columns refer to: the 
proper separation between the centers of the stars d/M ADM ; the baryon mass Mb of each star in solar masses; the total ADM mass M ADM 
in solar masses, as measured on the finite-difference grid; the approximate quasi-normal mode frequency of the fundamental mode / QNM in 
kHz; the SNR for Virgo, LIGO, Advanced LIGO and GEO. 



Model 


d/M ADM 


M b (M @ ) 


M ABM (M ) 


/qnm ( kHz ) 


SNR (Virgo) 


SNR (LIGO) 


SNR (Adv. LIGO) 


SNR (GEO) 


1.46-45-P 


14.3 


1.456 


2.681 


7.3 


1.92 


1.33 


12.54 


0.57 


1.46-45-IF 


14.3 


1.456 


2.681 




2.08 


1.45 


13.52 


0.62 


1.62-45-P 


12.2 


1.625 


2.982 


6.7 


2.15 


1.48 


13.29 


0.63 


1.62-45-IF 


12.2 


1.625 


2.982 


7.0 


2.29 


1.57 


14.42 


0.67 


1.62-60-P 


16.8 


1.625 


2.987 


6.5 


3.97 


3.15 


35.52 


1.48 



expression lfl3llfl32ll 

/ QNM a 3.23 (jj^£_) [l-0.63(l-a)°- 3 ] kHz, (38) 

this frequency is / QNM — 6.7 kHz for the BH produced by this 
binary (cf. Table inn. 

The PSD for the high-mass ideal-fluid binary, on the other 
hand, is more complex, with the inspiral peak at / w 0.75 kHz 
being accompanied by a number of other peaks, the most 
prominent having a similar amplitude at / sw 1.75 kHz and 
/ ss 3 kHz. These additional peaks (and also the smaller ones 
between the two) are obviously related to the post-merger 
phase at t > 5 ms and, in particular, to the dynamics of the 
HMNS formed after the merger and especially to the dynam- 
ics of the cores of the two NSs, which merge and bounce sev- 
eral times before the HMNS collapses to a BH, producing a 
small peak at / « 4 kHz. Also in this case even the funda- 
mental QNM has a frequency / QNM ~ 7.0 kHz (cf. Table ED 
and is therefore outside the range shown in Fig. [27] 

In a similar way it is possible to interpret the PSDs of the 
low-mass binaries. The polytropic one, in particular, shows 
an excess power at / r* 0.75 kHz due to the inspiral but also 
a very broad peak between / ~ 2 kHz and / w 3.5 kHz, that 
is related to the dynamics of the bar-deformed HMNS formed 
after the merger and persisting for several milliseconds. Also 
in this case a small excess power is seen at / > 4 kHz and is 
associated with the collapse to BH, whose fundamental QNM 
has a frequency / QNM — 7.3 kHz. Interestingly, the low- 
mass ideal-fluid PSD does not show the broad peak but a very 
narrow and high-amplitude one at / w 2kHz. This is obvi- 
ously related to the long-lived bar deformation of the HMNS, 
which we have followed for ~ 16 revolutions (as computed 
from the cycles of the quadrupolar gravitational radiation). 
At this stage it is unclear whether this prominent peak will 
survive when the simulations are repeated without the use 
of a 7r-symmetry and more conclusive results on this will 
be presented elsewhere ||90ll . Note that the high-frequency 
part of the PSD for the low-mass ideal-fluid binary (i.e. for 
/ > 2 kHz) is essentially zero, because of the absence of a 
collapse to BH, which for this binary takes place in an exces- 
sively long time. 

A fundamental piece of information necessary to assess 
the relevance of binary NSs as sources of gravitational waves 
comes from the calculation of the SNR which we have com- 
puted for interferometric detectors such as Virgo, LIGO, Ad- 
vanced LIGO and GEO. For all the models discussed above, 



including the high-mass polytropic binary with a larger initial 
separation of 60 km, the SNR has been computed as 

where Sh(f) is the noise power-spectral density for a given 
detector. The results computed assuming a source at a dis- 
tance of 10 Mpc are reported in Table [Till and show that, while 
a detection is ideally possible with the current interferometers 
[the SNR is 0(1)], it is unlikely in practice, given the small 
event rate at such distances, i.e. m O.Olyr -1 . On the other 
hand, larger SNRs O(10) can be obtained with advanced de- 
tectors. This also means that a detection of these sources up 
to a distance of 100 Mpc will be possible and so there will 
a higher event rate. Interestingly, binaries of the same mass, 
but described by a non-isentropic EOS have a slightly higher 
SNR and this is simply due to the increase in the delay for the 
collapse to BH. 

Both the small range in which the masses of NSs fall and 
the low sensitivity of present detectors in the high-frequency 
region, where a lot of the power is emitted, underline the im- 
portance of the inspiral phase for the detection. This is par- 
ticularly evident when comparing the large SNR of signals in 
which the inspiral is a significantly long part. The signal for 
the high-mass polytropic binary 1.62-60-P starts from an ini- 
tial separation of 60 km and spans over more than 5 orbits, 
resulting in a SNR which is a factor of 3 larger than the one of 
the other binaries, which have an initial separation of 45 km 
and merge in little more than 2 orbits. This result strongly mo- 
tivates the investigation, both through simulations and PN ap- 
proximations, of binaries inspiralling over timescales longer 
than the already long ones presented here. Stated differently, 
the study of longer simulations can be used to assess when the 
lower-order PN expressions are very accurate, while the study 
of the final part of the inspiral (say the last two orbits) can be 
used to determine those higher-order PN effects that have not 
been worked out analytically yet. 



V. CONCLUSIONS 

We have discussed accurate general-relativistic simulations 
of binary systems of equal-mass NSs which inspiral starting 
from irrotational configurations in quasi-circular orbit. Span- 
ning over ~ 30 ms, our simulations are the longest of their 
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kind and provide the first complete (within an idealized treat- 
ment of the matter) description of the inspiral and merger of 
a NS binary leading to the prompt or delayed formation of a 
BH and to its ringdown. 

More specifically, we have considered binary NSs with two 
different initial masses: low-mass binaries with M ADM = 
2.681M Q and high-mass binaries with M ADM = 2.982M©. 
Such binaries have then been evolved using two different 
EOSs: namely an isentropic {i.e. polytropic) EOS and a non- 
isentropic (i.e. ideal-fluid) EOS. Despite the use of only sim- 
ple, analytical EOSs, we were able to reproduce some of the 
aspects that a more realistic EOS would yield. In particular, 
we have shown that the polytropic EOS leads either to the 
prompt formation of a rapidly rotating BH surrounded by a 
dense torus in the high-mass case, or, in the low-mass case, to 
a HMNS which develops a bar, emits large amounts of gravita- 
tional radiation and eventually experiences a delayed collapse 
to BH. Conversely, we have shown that the ideal-fluid EOS 
inevitably leads to a further delay in the collapse to BH as a 
result of the larger pressure support provided by the temper- 
ature increase via shocks. In this case the temperature in the 
formed HMNS can reach values as high as 10 11 — 10 12 K, so 
that the subsequent dynamics and especially the time of the 
collapse can be reduced if cooling mechanisms, such as the 
direct-URCA process, take place. 

With the exception of the low-mass ideal-fluid binary, 
whose HMNS is expected to collapse to BH on a timescale 
which is computationally prohibitive (i.e. ~ 110 ms), all the 
binaries considered lead to the formation of a BH surrounded 
by a rapidly rotating torus. The masses and dimensions of the 
tori depend on the EOS, but are generically larger than those 
reported in previous independent studies, with masses up to 
w O.O7M . Confirming what was reported in ref. JIU], we 
have found that the amount of angular momentum lost during 
the inspiral phase can influence the mass of the torus for bi- 
naries that have the same EOS. In particular, the models that 
lose less angular momentum during the inspiral, the compar- 
atively low-mass binaries, are expected to have comparatively 
high-mass tori. A more detailed study of the dynamics of the 
torus (especially when produced from non-equal-mass bina- 
ries) and of its implication for short hard GRBs will be the 
subject of a following paper l90ll . 

Most of the binaries considered have an initial coordinate 
separation of 45 km and merge after ~ 2 orbits or, equiva- 
lent^, after ~ 6 ms. However, we have also considered a 
high-mass polytropic binary with an initial coordinate sepa- 
ration of 60 km, which merges after ~ 5 orbits or, equiva- 
lent^, after ~ 20 ms. As a stringent test of the accuracy of 
our results we have carried out a systematic comparison be- 
tween identical binaries starting at different initial separations. 
Such a comparison, which has never been performed before, 
has shown that there is an excellent agreement in the inspiral 
phase (as expected from the lowest-order PN approximations), 
but also small differences at the merger and in the subsequent 
evolution. These results provide us with confidence on our 
ability to perform long-term accurate simulations of the inspi- 
ral phase, and also open the prospect of investigating higher- 
order PN corrections. 



Besides the study of the bulk dynamics of the two NSs, 
we have also investigated the small-scale hydrodynamics of 
the merger and the possibility that dynamical instabilities 
develop. In this way we have provided the first quantita- 
tive description of the onset and development of the Kelvin- 
Helmholtz instability, which takes place during the first stages 
of the merger phase, when the outer layers of the stars come 
into contact and a shear interface forms. The instability curls 
the interface forming a series of vortices which we were able 
to resolve accurately using the higher resolutions provided by 
the AMR techniques. Since the development of this instability 
is essentially independent of the EOS used or of the masses of 
the NSs, it could have important consequences in the genera- 
tion of large magnetic fields. Also this aspect will be further 
investigated in a subsequent work Ir90ll . 

Given the importance of binary NSs as sources of gravita- 
tional waves, special attention in this work has been dedicated 
to the analysis of the waveforms produced and to their prop- 
erties for the different configurations. In particular, we have 
found that the largest loss rates of energy and angular momen- 
tum via gravitational radiation develop at the time of the col- 
lapse to BH and during the first stages of the subsequent ring- 
down. Nevertheless, the configurations which emit the high- 
est amount of energy and angular momentum are those with 
lower masses, since they do not collapse promptly to a BH. In- 
stead they produce a violently oscillating HMNS, which emits 
copious gravitational radiation, while rearranging its angular- 
momentum distribution, until the advent of the collapse to BH. 
We have also found that although the gravitational-wave emis- 
sion from NS binaries has spectral distributions with large 
powers at high frequencies (i.e. f > 1 kHz), a signal-to-noise 
ratio as large as 3 can be estimated for a source at lOMpc 
using the sensitivity of currently operating gravitational-wave 
interferometric detectors. 

Several aspects of the simulations reported here could be 
improved and probably the most urgent among them are the 
use of more realistic EOSs and the inclusion of magnetic 
fields via the solution of the MHD equations. Recent calcu- 
lations l48l l49ll have indeed shown that the corrections pro- 
duced by strong magnetic fields could be large and are proba- 
bly very likely to be present. Work is in progress towards these 
improvements using the code developed in ref. f5lll . The re- 
sults of these investigations will be presented in forthcoming 
works. 



APPENDIX A: CHARACTERIZING THE TRUNCATION 
ERROR 

1. The influence of numerical methods 

The inherent numerical viscosity of the numerical method 
used for the reconstruction of the variables on cell interfaces 
is crucial to determine the time of the merger. As one might 
expect, lower-order reconstruction schemes result in an an- 
ticipated merger due to their higher numerical viscosity, as 
Fig. [28] shows (for a review of the numerical methods imple- 
mented in Whisky, see III CI ). The results in the test Simula- 
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FIG. 28: Comparison of the rest-mass density normalised to its value 
at t = for evolutions performed with different numerical methods; 
the solid line refers to an evolution performed using the Marquina 
flux formula and a PPM reconstruction, the dotted line to HLLE- 
PPM and the dashed line to HLLE-TVD (van Leer slope limiter). 
These data refer to an initial configuration not present in Table|l](see 
text for details) and to an evolution with the ideal-fluid EOS. 



tions presented in the figure were produced through the evo- 
lution of initial data that are not listed in Table U i.e. proper 
separation between the centers of the stars d/M AT)M = 12.6; 
baryon mass of each star Mb = 1.78 M©; total ADM mass 
A/ adm = 3.24 M ; angular momentum J = 9.93 Afg = 
8.75 x 10 49 gcm 2 /s; initial orbital angular velocity Hq = 
9.39 x 1CP 3 = 1.9rad/ms; approximate mean radius of each 
star R = 8.4 Af = 12 km; ratio of the polar to the equatorial 
coordinate radius of each star r p /r e = 0.945. 

In particular, in Fig. [28] we show the differences in the evo- 
lution of the rest-mass density normalized to its initial value 
when different numerical methods are used for the evolution: 
the solid line refers to an evolution performed using the Mar- 
quina flux formula and a PPM reconstruction (which is our 
usual choice), the dotted line to the HLLE approximate Rie- 
mann solver with PPM reconstruction and the dashed line to 
the HLLE solver with TVD reconstruction (in particular, the 
van Leer slope limiter was used). Smaller changes in the 
merger time and in the evolution of the HMNS are observed 
also by changing some parameters of the PPM reconstruction 
method, in particular those related to the shock detection, that 
is the parameters that define how big a jump in the evolved 
variable has to be, in order to be considered a discontinuity 
and treated as such. 

We have found instead that the choice of approximate Rie- 
mann solver does not influence significantly the evolution of 
the coalescence. As one can see from Fig. |28j when cou- 
pled with the PPM reconstruction, both the Marquina and the 
HLLE solvers produce very similar dynamics and the time of 



FIG. 29: Comparison of the time evolution of the coordinate separa- 
tion (upper panel) and the proper separation (lower panel) between 
the stellar centres in case the initial Meudon shift is used (dashed 
line) and in case the initial shift is set to zero (continuous line). 

the merger is almost the same. The situation changes when a 
lower-order reconstruction method, such as the van Leer one, 
is used. In this case the numerical viscosity is large and the 
time of the merger is very different, i.e. ~ 4 ms instead of 
ss 6.5 ms. 

From these tests one can then learn that the numerical vis- 
cosity of the evolution method is very important in this sce- 
nario, being responsible for changes in the dynamics and also 
in the estimate of the gravitational-wave emission. Of course, 
one should always employ the least viscous method available. 



2. The influence of the initial gauge conditions 

We have found that using the shift profile given in the 
Meudon data introduces a considerable amount of gauge dy- 
namics, which can be avoided by setting the initial shift to 
zero. We recall that the Meudon shift condition is determined 
through the Killing equation which is implicit in the quasi- 
equilibrium assumption for binary systems ll2uTl . A clear way 
to highlight this feature is a comparison of the time evolution 
of the coordinate separation between the stellar centres. This 
is shown in Fig. |29j which offers a comparison of the time 
evolution of the coordinate separation (upper panel) and the 
proper separation (lower panel) between the stellar centres in 
case the initial Meudon shift is used (dashed line) and in case 
the initial shift is set to zero (continuous line). The evolution 
equation for the shift is the same for the two simulations. 

It is clear that the coordinate orbit of the evolution started 
with the Meudon shift has a noticeable amount of eccentricity 
(which appears as large oscillations of the coordinate separa- 
tion of the stars during the inspiral), which is absent in the 
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simulation in which the shift is zero at the initial time. In ad- 
dition, the adoption of a zero initial shift results into a larger 
initial violation of the 2-norm of the Hamiltonian constraint, 
which is however below 10~ 6 for the typical resolution used. 

As a final remark we note that the proper separations of the 
stars, the maximum of the rest-mass density and other gauge- 
invariant quantities, such as the gravitational waveforms, are 
instead very similar during the inspiral phase. 
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